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Abstract: 

We  study  the  stability  properties  of,  and  the  phase  error  present  in,  a  finite  element  scheme 
for  Maxwell’s  equations  coupled  with  Debye  or  Lorentz  polarization.  In  one  dimension 
we  consider  a  second  order  formulation  for  the  electric  field  with  an  ordinary  differential 
equation  for  the  polarization  added  as  an  auxiliary  constraint.  The  finite  element  method 
uses  linear  finite  elements  in  space  for  the  electric  field  as  well  as  the  polarization,  and 
a  theta  scheme  for  the  time  discretization.  Numerical  experiments  suggest  the  method  is 
unconditionally  stable  for  both  Debye  and  Lorentz  models.  We  compare  the  stability  and 
phase  error  properties  of  the  method  presented  here  with  those  of  finite  difference  methods 
that  have  been  analyzed  in  the  literature.  We  also  conduct  numerical  simulations  that  verify 
the  stability  and  dispersion  properties  of  the  scheme. 
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1  Introduction 


Noninvasive  interrogation  of  the  interior  of  tissues  and  other  materials  by  electromagnetic 
waves  has  important  applications  in  various  fields  including  medical  imaging  for  the  early 
detection  of  anomalies  and  nondestructive  damage  detection  in  aircraft  [2,  3].  For  example, 
microwave  imaging  for  breast  cancer  detection  is  expected  to  be  safe  for  the  patient  and 
has  the  potential  to  detect  very  small  cancerous  tumors  in  the  breast  [16].  This  ability  for 
detection  is  based  on  the  difference  in  electrical  properties  of  malignant  and  normal  tissues. 
Biological  tissue  interactions  with  the  fields  are  defined  by  their  complex  permittivity  which 
is  a  function  of  the  various  electric  and  magnetic  polarization  mechanisms  and  conductivity 
of  the  biological  medium.  Similarly,  nondestructive  damage  detection  in  materials  for  the 
detection  of  defects  such  as  cracks  is  based  upon  the  changes  in  the  electrical  properties 
that  occur  due  to  the  presence  of  these  defects.  Thus,  one  of  the  aims  of  electromagnetic 
interrogation  problems  is  the  determination  of  the  dielectric  properties  of  the  materials 
under  investigation. 

To  computationally  simulate  these  electromagnetic  interrogation  problems  requires  set¬ 
ting  up  a  suitable  inverse  problem  which  involves  numerous  forward  simulations  of  the 
propagation  of  transient  electromagetic  waves  in  lossy  dispersive  dielectrics  such  as  bio¬ 
logical  tissue.  Hence,  the  development  of  accurate,  consistent  and  stable  discrete  forward 
solvers  is  very  important,  as  errors  in  the  numerical  solvers  can  result  in  inaccurate  deter¬ 
mination  of  the  dielectric  properties  that  determine  the  characteristics  of  the  material  being 
investigated. 

The  electric  and  magnetic  fields  inside  a  material  are  governed  by  the  macroscopic 
Maxwell’s  equations  along  with  constitutive  laws  that  account  for  the  response  of  the  ma¬ 
terial  to  the  electromagnetic  field.  In  special  cases  Maxwell’s  equations  can  be  reduced  to  a 
vector  wave  equation  in  the  electric  or  magnetic  fields.  Numerical  approximation  algorithms 
of  time- dependent  wave  equations  and  Maxwell’s  equations  introduce  error  into  the  ampli¬ 
tude  and  speed  of  the  propagating  waves.  These  errors  include  dissipation ,  the  dampening 
of  some  frequency  modes,  and  dispersion ,  the  frequency  dependence  of  the  phase  velocity 
of  numerical  wave  modes  in  the  computational  grid. 

Dielectric  materials  have  actual  physical  dispersion.  The  complex  electric  permittivity  of 
a  dielectric  medium  is  frequency  dependent  (has  dielectric  dispersion).  Thus,  an  appropriate 
discretization  method  should  have  a  numerical  dispersion  that  matches  the  model  dispersion 
as  closely  as  possible.  Dielectric  materials  also  have  physical  dissipation,  or  attenuation, 
which  must  also  be  correctly  computed  by  a  numerical  method.  In  particular,  if  a  method 
does  not  sufficiently  damp  initial  disturbances,  possibly  due  to  round-off  error,  the  method 
can  become  unstable.  Certain  algorithms  have  criteria  based  on  discretization  parameters 
to  determine  when  it  may  be  unstable;  these  are  called  conditionally  stable  methods. 

The  stability  and  dispersion  properties  for  the  finite  difference  time  domain  (FDTD) 
schemes  applied  to  Maxwell’s  equations  in  free  space  are  well  known  (see  [35]).  Addition¬ 
ally,  different  time  domain  finite  element  methods  have  also  been  devised  for  the  numerical 
approximation  of  Maxwell’s  equations  in  free  space  (see  [27,  23]  and  the  references  therein). 
While  free  space  analyses  for  finite  element  and  finite  difference  methods  are  well  doc¬ 
umented,  stability  and  phase  error  analysis  for  dispersive  dielectrics  has  been  primarily 
focused  on  finite  difference  methods  (see  [28,  35]  for  standard  second  order  methods  and 
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[29]  for  higher  order  schemes).  The  treatment  for  finite  element  methods  has  been  limited  to 
scalar-potential  formulations  to  model  dielectric  dispersion  at  low  frequencies  ([31]),  scalar 
Helmholtz  equation  ([24]),  and  in  some  cases,  hybrid  methods  ([15]). 

There  are  several  FDTD  extensions  that  have  been  developed  to  model  electromagnetic 
pulse  propagation  in  dispersive  media.  One  technique  is  to  add  to  Maxwell’s  equations  a 
set  of  ODE’s  that  relate  the  electric  displacement  D{t)  to  the  electric  field  E(t)  [18],  or 
a  set  of  ODE’s  that  model  the  dynamic  evolution  of  the  polarization  vector  P(t)  driven 
by  the  electric  field  [20,  19].  Dielectric  dispersion  can  be  expressed  in  the  time  domain 
as  a  convolution  integral  involving  the  electric  field  and  a  causal  susceptibility  function. 
The  recursive  convolution  method  [26,  25,  21]  uses  a  recursive  technique  to  update  the 
convolution  representation  of  the  constitutive  law  along  with  the  FDTD  time  update  of 
Maxwell’s  equations.  There  are  other  methods  such  as  the  Z-transform  [34,  33]  and  the 
TLM  method  [9]  that  have  also  been  used  to  model  pulse  propagation  in  dispersive  media. 
Many  of  these  methods  have  been  compared  and  analyzed  for  their  numerical  errors  and 
stability  properties  [28,  36,  30,  10,  12]. 

In  this  paper  we  study  the  propagation  characteristics  of  the  discretized  Maxwell’s  equa¬ 
tions  with  Debye  or  Lorentz  polarization,  using  a  finite  element  method,  in  terms  of  numer¬ 
ical  stability  and  dispersion  analyses.  We  obtain  information  about  the  expected  accuracy 
of  the  method  from  the  construction  of  the  dispersion  relation  which  relates  the  numerical 
wave  number  k  to  the  frequency  u  for  waves  propagating  in  the  finite  element  grid  and  then 
compare  with  the  dispersion  relation  for  the  corresponding  continuous  model  (differential 
equations).  We  compare  the  stability  and  dispersion  properties  of  the  finite  element  method 
with  those  of  finite  difference  methods  analyzed  in  [28].  Such  finite  element  methods  have 
been  used  for  the  electromagnetic  interrogation  of  dielectric  media  with  a  metal  backing, 
for  the  determination  of  dielectric  properties  and  geometrical  dimensions  of  the  medium 
[5],  as  well  as  for  the  detection  of  cracks  in  composite  materials  [8].  In  [1],  finite  element 
methods  were  used  to  interrogate  dielectric  media  with  acoustic  waves  as  virtual  reflectors. 

We  analyze  one  dimensional  models  to  which  we  apply  standard  linear  finite  elements 
for  the  spatial  discretization  of  the  electric  field  and  the  polarization.  Maxwell’s  equations 
are  reduced  to  a  wave  equation  in  the  electric  field,  and  the  constitutive  law  in  the  medium 
involves  an  ordinary  differential  equation  that  describes  the  dynamic  evolution  of  the  po¬ 
larization  driven  by  the  electric  field.  The  entire  system  is  rewritten  in  first  order  form  and 
time  discretized  using  a  theta  scheme. 

The  numerical  stability  results  for  Debye  and  Lorentz  media  suggest  that  the  extension 
of  the  finite  element  scheme  for  these  media  retain  the  unconditional  stability  property  of 
the  scheme  in  free  space.  The  well  known  Lax-Richtmayer  theorem  [32]  states  that  the 
convergence  of  consistent  difference  schemes  to  initial  value  problems  represented  by  PDE’s 
is  equivalent  to  stability.  Thus,  unconditional  stability  is  a  desirable  property  of  difference 
schemes  as  it  allows  the  choice  of  the  time  step  to  be  determined  by  the  physical  dimensions 
of  the  problem,  such  as  a  relaxation  time.  The  FDTD  methods,  on  the  other  hand,  are 
conditionally  stable.  Since  the  stability  condition  is  determined  by  the  smallest  cell  size  in 
the  domain,  the  FDTD  analysis  of  very  fine  geometric  structures  requires  a  large  number 
of  time  iterations.  In  the  finite  element  case  unconditional  stability  allows  us  to  choose 
the  Courant  number ,  which  relates  the  time  step  to  the  mesh  step  size,  to  minimize  the 
numerical  phase  error. 
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The  numerical  dispersion  analysis  shows  that  for  accuracy  we  need  to  resolve  the  shortest 
time  scale  in  the  problem  which  agrees  with  the  result  obtained  in  [28].  This  is  reflected 
in  the  fact  that  the  time  step  for  simulating  Debye  media  should  be  chosen  to  be  about 
O(10~3)t,  where  r  is  the  relaxation  time  of  the  medium,  whereas  for  Lorentz  media,  the 
time  step  should  be  chosen  to  be  the  minimum  of  O(10~2)t  and  0(10~2)(^),  where  loq  is 
the  resonance  frequency  of  the  medium. 

We  thus  determine  a  guideline  for  users  in  which  the  time  step  is  chosen  to  minimize 
the  dissipation  of  the  scheme  and  the  Courant  factor  is  chosen  to  minimize  the  phase  error, 
providing  good  agreement  between  the  exact  and  numerical  complex  permittivity. 

2  Model  formulation 

We  consider  Maxwell’s  equations  which  govern  the  electric  field  E  and  the  magnetic  field 
H  in  a  domain  O  from  time  0  to  T  given  as 


(J I  J  — *  — »  _ 

+  J  —  V  x  H  =  0  in  (0,  T)  x  O,  (2.1a) 

dB  -> 

-7^  +  V  x  E  =  0  in  (0,  T)  x  O,  (2-lb) 

V  •  D  =  0  in  (0,T)  x  O,  (2.1c) 

V  •  B  =  0  in  (0,T)  x  SI,  (2.1d) 

E(0,x)  =  0  in  O,  (2.1e) 

H(0,x)  =  0  in  O.  (2. If ) 


The  fields  D,  B  are  the  electric  and  magnetic  flux  densities  respectively.  All  the  fields  in 
(2.1)  are  functions  of  position  x  =  (x,  y,  z )  and  time  t.  We  have  J  =  Jc  +  Js,  where  Jc  is  a 
conduction  current  density  and  Js  is  the  source  current  density.  However,  we  will  assume 
Jc  =  0  in  this  paper,  as  we  are  interested  in  dielectrics  with  no  free  charges.  Appropriate 
boundary  conditions  are  added  to  system  (2.1)  to  terminate  the  computational  domain. 

Constitutive  relations  which  relate  the  electric  and  magnetic  fluxes  D,  B  to  the  electric 
and  magnetic  fields  E,  H  are  added  to  these  equations  to  make  the  system  fully  determined 
and  to  describe  the  response  of  a  material  to  the  electromagnetic  fields.  In  free  space,  these 
constitutive  relations  are  D  =  eo E,  and  B  =  where  eo  and  /xq  are  the  permittivity 

and  the  permeability  of  free  space,  respectively,  and  are  constant  [17].  In  general  there  are 
different  possible  forms  for  these  constitutive  relationships.  In  a  frequency  domain  formula¬ 
tion  of  Maxwell’s  equations,  these  are  usually  converted  to  linear  relationships  between  the 
dependent  and  independent  quantities  with  frequency  dependent  coefficient  parameters.  We 
will  consider  the  case  of  a  dielectric  in  which  magnetic  effects  are  negligible.  Thus,  within 
the  dielectric  medium  we  have  constitutive  relations  that  relate  the  flux  densities  D,  B  to 
the  electric  and  magnetic  fields,  respectively,  as 


D  =  e0E  +  P. 
B  =  n0H. 


4 


(2.2a) 

(2.2b) 


In  (2.2a),  the  quantity  P  is  called  the  macroscopic  electric  polarization.  (A  discussion  of  the 
relationship  between  the  macroscopic  polarization  and  the  microscopic  material  properties 
leading  to  distributions  of  relaxation  times  and  other  dielectric  parameters  in  the  constiu- 
tive  laws  can  be  found  in  [6].)  Electric  polarization  may  be  defined  as  the  electric  field 
induced  disturbance  of  the  charge  distribution  in  a  region.  This  polarization  may  have  an 
instantaneous  component  as  well  as  delayed  effects;  the  latter  will  usually  have  associated 
time  constants  called  relaxation  times  which  are  denoted  by  r.  We  define  the  instantaneous 
component  of  the  polarization  to  be  related  to  the  electric  field  by  means  of  a  dielectric 
constant  eo  and  a  susceptibility  y.  The  remainder  of  the  electric  polarization,  called  the 
relaxation  polarization,  is  denoted  as  Pr.  Therefore,  we  have 

P  =  Pi  +  Pr  =  (oX.P  +  Pr,  (2-3) 

and  hence  the  constitutive  law  (2.2a)  becomes 

D  =  eo  erE  +  Pr,  (2.4) 

where  er  =  (1  +  x)  is  the  relative  permittivity  of  the  dielectric  medium.  We  will  henceforth 
denote  Pr  by  P,  as  the  instantaneous  polarization  will  be  absorbed  into  the  dielectric 
constant  er.  The  following  section  defines  the  equations  for  polarization  models  of  interest 
in  this  paper. 

2.1  Models  for  Polarization 

To  describe  the  behavior  of  the  media’s  relaxation  polarization  P,  one  may  use  ordinary 
differential  equation  models  that  describe  microscopic  polarization  mechanisms  such  as 
dipole  or  orientational  polarization  (Debye),  as  well  as  ionic  and  electronic  polarization 
(Lorentz),  and  other  frequency  dependent  polarization  mechanisms  [4].  For  more  complex 
dielectric  materials,  a  simple  Debye  or  Lorentz  polarization  model  is  often  not  adequate  to 
characterize  the  dispersive  behaviour  of  the  material.  One  can  then  turn  to  combinations 
of  Debye,  Lorentz,  or  even  more  general  nth  order  mechanisms  [5]  as  well  as  Cole-Cole 
type  (fractional  order  derivative)  models  [11].  Additionally,  materials  may  be  represented 
by  a  distribution  of  the  associated  time  constants  or  even  a  distribution  of  polarization 
mechanisms  (see  [7,  6]).  In  this  report  we  concentrate  our  analysis  on  single  pole  Debye 
and  Lorentz  polarization  models. 

2.1.1  Debye  Model 

The  differential  equation  of  the  Debye  model  for  orientational  or  dipolar  polarization  is  given 
by 

dP 

t— +  P  =  e0(es  -  e^E.  (2.5) 

Here  es  is  the  static  relative  permittivity.  The  presence  of  instantaneous  polarization  is 
accounted  for  in  this  case  by  the  coefficient  in  the  electric  flux  equation.  That  is, 
er  =  £oo-  The  remainder  of  the  electric  polarization  is  seen  to  be  a  decaying  exponential 
with  relaxation  parameter  r,  driven  by  the  electric  field,  less  the  part  included  in  the 
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instantaneous  polarization.  This  model  was  first  proposed  by  Debye  [14]  to  model  the 
behaviour  of  materials  that  possess  permanent  dipole  moments.  The  magnitude  of  the 
polarization  term  P  represents  the  degree  of  alignment  of  these  individual  moments  and  is 
based  on  a  uniformity  assumption  at  the  molecular  level  (see  [6] ) .  The  choice  of  coefficients 
in  (2.5)  gives  a  physical  interpretation  to  es  and  as  the  relative  permittivities  of  the 
medium  in  the  limit  of  the  static  field  and  very  high  frequencies,  respectively.  In  the  static 
case,  we  have  Pt  =  0,  so  that  P  =  eo(es  —  Coo )E  and  D  =  eo esE.  For  very  high  frequencies, 
rPt  dominates  P  so  that  P  «  0  and  D  =  e^e^E  (thus  the  notation  of  oo). 

The  Debye  model  is  most  often  used  to  model  electromagnetic  wave  interactions  with 
water-based  substances,  such  as  biological  materials.  In  particular,  biological  tissue  is  well 
represented  by  multi-pole  Debye  models,  by  accounting  for  permanent  dipole  moments  in 
the  water.  The  Debye  model  has  other  physical  characteristics  which  make  it  attractive 
from  an  analytical  point  of  view  (for  details,  see  [36]). 


2.1.2  Lorentz  Model 

The  Lorentz  model  for  electronic  polarization  in  differential  form  is  represented  with  the 
second  order  equation: 

o2p  i  f)p 

-dP+r~di+^P  =  e^E'  (2'6) 

where  ed  =  es  —  e ^  and  loq  is  the  resonance  frequency  of  the  material. 

The  Lorentz  model  is  formulated  by  modeling  the  atomic  structure  of  the  material  as  a 
damped  vibrating  system.  Applying  classical  Newtonian  laws  of  motion,  the  displacement 
of  the  outermost  shell  of  the  atom  is  found  to  be  a  solution  to  a  second-order,  ordinary 
differential  equation  [36]. 


2.2  Reduction  to  One  Dimension 


The  electric  field  is  assumed  to  be  polarized  to  have  oscillations  in  the  x-z  plane  only,  as 
described  in  [5].  Restricting  the  problem  to  one  dimension,  we  can  write  the  electric  and 
magnetic  fields,  E  and  H  respectively,  as  follows 

E(t,  x)  =  iE(t,  z) 

H(t,x)  =  jH(t,z), 


so  that  we  are  only  concerned  with  the  scalar  values  E(t,z )  and  H(t,z).  In  this  case 
Maxwell’s  equations  become: 


dE  dH 

~  ~^°~dt 
dH  _dD 

~~dt+  s 


(2.7a) 

(2.7b) 


We  take  the  partial  derivative  of  Equation  (2.7a)  with  respect  to  z,  and  the  partial  of 

a2  t  t 

Equation  (2.7b)  with  respect  to  t.  Equating  the  terms  in  each,  and  thus  eliminating 
the  magnetic  field  H,  we  have: 


E" 


Ho  [D  +  J& 
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(where  '  denotes  2  derivatives  and  ’  denotes  time  derivatives).  Using  the  constitutive  law 
for  he  electric  flux  density  given  by  D  =  eoe^E  +  P,  we  have 

WUoeoo E  +  P  -  E"  =  -/i0  js  in  U  =  [a,  b] .  (2.8) 


In  order  to  have  a  finite  computational  domain,  we  impose  absorbing  boundary  conditions 
at  z  =  a  and  z  =  b,  which  are  modeled  as 


E-cE1 


=  0; 


E  +  cE' 


-  z=b 


=  0, 


(2.9) 


where  c  =  1/y/eopo  is  the  speed  of  light  in  vacuum.  With  these  boundary  conditions, 
any  incident  signal  passes  out  of  the  computational  domain,  and  does  not  return.  The 
homogeneous  initial  conditions  in  ID  become 


E( 0,  z)  =  0,  P( 0,  z)  =  0,  P(0,  z)  =  0,  E( 0,  z)  =  0. 


(2.10) 


3  Numerical  Solution 

3.1  Spatial  Discretization  Using  Finite  Elements 

We  apply  a  finite  element  method  using  standard  piecewise  linear  one  dimensional  basis 
elements  to  discretize  the  model  (2.8)  in  space.  Let  N  be  the  number  of  intervals  in  the 
discretization  of  z,  and  A2:  =  (6  —  a)/N,  then  the  finite  element  discretization  has  an  order 
of  accuracy  of  0(Az2).  The  resulting  system  of  ordinary  differential  equations  after  the 
spatial  discretization  is  the  semi-discrete  form 


with  either 


or 


eocHoeoMe  +  /jlq  Mp  +  Be  +  Ke  =  no  J, 

p  +  Xp  =  eofirfAe,  (Debye  Media), 
p  +  Xp  +  loqP  =  eo e,  (Lorentz  Media), 


(3.1) 

(3.2) 

(3.3) 


where  A  =  1.  The  vectors  e  and  p  represent  the  values  of  E,  and  P,  respectively  at  the 
nodes  z*  =  iAz.  The  mass  matrix  M  has  entries 

Mij  =  (0i,  4>j)  :=  /  ( h<t>jdz , 


where  are  the  basis  functions.  The  stability  matrix  K  has  entries 

,,  -  f 


Kij  =  (<l>i,<t>'j)  ■=  I  4>i4>'jdz. 


The  matrix  B  results  from  the  boundary  conditions  where 

Bij  =  ^  [<f>i(a)<j)j(a)  -  ^(6)^(6)] . 
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Finally,  J  is  defined  as 


Ji  —  ( ■  Ts )  • —  I  Js4>idz. 

J  a 

For  Debye  media,  we  differentiate  (3.2)  and  substitute  for  p  into  (3.1)  to  obtain  an 
equation  only  dependent  explicitly  on  p,  given  as 


A2 


1 


£oo Me  +  (-B  +  €d\M)e  +  (c  K  —  edX  M)e  H - Mp  —  — J , 

eo  eo 

V  +  \p  -  e0erfAe  =  0. 

In  the  above  c  =  l/^/eo^o  is  the  speed  of  light  in  vacuum. 

For  Lorentz  media  we  substitute  (3.3)  into  (3.1)  to  get 


Wn 


A 


1 


(3.4a) 

(3.4b) 


too  Me  +  Be  +  ( czK  +  qwjM)e - -  Mp - Mp  =  — J, 

eo  eo  eo 

p  +  Xp  +  wlp  =  eoedu>le. 

For  linear  finite  elements  in  one  dimension,  the  entries  of  the  mass  matrix  M  are 


(3.5a) 

(3.5b) 


Mij  = 


2Az/3,  if  0  <  i  =  j  <  N, 
Az/3,  else  if  i  =  j  =  0  or  N, 
k  Az/6,  else  if  i  =  j  ±  1. 


(3.6) 


The  entries  of  the  stiffness  matrix  can  be  calculated  as 


Kij  —  ' 


2/Az,  if  0  <  i  =  j  <  N, 

1/A z,  else  if  *  =  j  =  0  or  N, 
—  1/Az,  else  if  i  =  j  ±  1. 


(3.7) 


3.2  Time  Discretization  Using  Finite  Differences 

In  order  to  solve  the  semi-discrete  form  of  our  equations  we  may  convert  each  coupled  second 
order  system  of  equations  into  one  larger  first  order  system  and  apply  a  theta  method  [22] . 
As  our  finite  element  method  is  second  order  in  space,  if  we  choose  6  =  ^  the  discretization 
is  second  order  in  time  as  well  (thus,  we  have  used  6  =  \  throughout).  Therefore,  for 
appropriately  smooth  data  and  with  At  =  O(A^),  the  combined  method  is  second  order  in 
time  and  space.  For  the  system  of  equations  corresponding  to  the  Debye  model  we  consider 
an  additional  method  of  temporal  finite  difference  for  comparison. 


3.2.1  Debye  Media 

In  the  first  method  to  be  considered,  we  will  convert  the  coupled  second  order  system  of 
equations  into  one  larger  first  order  system  and  simply  apply  a  theta  method  as  described 
above.  In  the  second  method,  we  will  solve  first  for  the  polarization  with  a  forward  differenc¬ 
ing  scheme  using  the  initial  conditions,  and  then  use  polarization  solution  vector  to  update 


a  second  order  central  difference  scheme  for  the  magnitude  of  the  electric  field.  We  then 
continue  this  process  iteratively,  alternating  between  solving  for  p  and  for  e.  Both  methods 
are  second  order  in  time  and  space  for  appropriately  smooth  data  (and  with  At  =  O(h)). 


Method  1:  For  Debye  media  we  convert  (3.4a)-(3.4b)  into  a  first  order  system  of  equations 
in  three  unknowns,  X  =  [e,p,  d]T,  where  d  =  e.  resulting  in 


with 


and 


I<  = 


MX  +  KX  =  J, 


M  = 

0 

— €o 


I  0  0 

0  I  0 

0  0  eoo  M 
0 


A 


-I 

0 


(c2 K  -  ed\2M)  B  +  ed\NI 


J=[0  0  ij]  . 

We  apply  a  theta-schenre  to  (3.8)  to  obtain 

(M  +  9AtK)Xn+1  =  (M  -  (1  -  9)AtK)Xn  +  ( 0Jn+l  +  (1  -  9)  Jn). 


(3.8) 


(3.9) 


(3.10) 

(3.11) 

(3.12) 


For  9  =  0.5  the  scheme  can  be  written  as 

(M  +  ^ K)Xn+1  =  (M  -  ^ K)Xn  +  \{Jn+l  +  Jn)-  (3.13) 

where,  Jn  =  J(nAt)  and  Xn  =  X(nAt). 

As  we  are  assuming  a  fixed  time  step  At,  the  matrix  to  be  inverted  does  not  change 
over  time.  Therefore,  we  solve  for  the  LU  factorization  at  the  beginning  and  use  back 
substitution  at  each  time  step.  Since  the  matrices  M  and  K  are  sparse,  we  employ 
the  sparse  matrix  package  UMFPACK  [13].  For  problems  where  the  number  of  time 
steps  is  small  compared  with  the  number  of  nodes,  it  may  be  more  efficient  to  use  an 
iterative  solver  with  the  previous  solution  vector  as  the  initial  iterate. 

Method  2:  In  our  second  method  we  use  a  second  order  central  difference  scheme  to  solve 
(3.4a).  Our  approach  is  to  first  solve  for  p  using  a  0-method,  and  then  use  that  approx¬ 
imation  to  solve  for  e  at  the  next  time  step.  Thus,  our  finite  difference  approximation 
for  (3.4b)  is 

pn+ 1  =pn-  At\pn+e  +  At\eden+e  (3.14) 

where  vn+e  =  9vn+1  +  (1  —  9)vn,  for  v  =  e  or  v  =  p.  This  implies 

p”+1=!’”  +  r TJKFe^-^-  (3J5) 
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Once  we  have  pn+l  we  can  solve  for  en+2.  Applying  second  order  central  difference 
with  averaging  to  (3.4a)  gives 


^eoo M(en+2  -  2en+1  +  en )  +  +  edXM)(en+2  -  en) 

+  ~(c2K  -  edX2M)(en+2  +  2en+1  +  en )  =  -Jn+1  -  —Mpn+1. 
4  eo  eo 


(3.16) 


Defining  hT  =  AA t  and  solving  for  the  en+2  term  we  have 


*  _  M  +  1a tB 

2  4/2 


en+2  = 


edh2 , 


(2eoo  H - Y~)M  ~  2 


At2 


-K 


n+ 1 


At2  /?2 

en  +  jn+l  _  -^Mpn+1 

eo  £0 


(3.17) 


or  equivalently, 

At2  /)2 

Aie^2  =  A2en+1  +  A3en  +  —  Jn+1  -  -^Mpn+1.  (3.18) 

eo  e0 

Note  that  in  this  case  A\  is  tridiagonal  and  the  matrix  is  the  same  for  each  time  step, 
so  we  may  store  the  Crout  LU  factorization  and  use  back  substitution  to  solve  the 
system  at  each  time  step.  For  tridiagonal  matrices  the  factorization  and  the  back 
substitution  are  both  order  O(N)  [22], 

Again,  for  6  =  |,  (3.15)  will  be  second  order  in  time  if  the  corresponding  solution  is 
C3  in  time.  Equation  (3.18)  is  also  second  order  in  time  assuming  an  exact  solution 
for  P,  and  that  E  has  four  continuous  time  derivatives  (for  the  second  order  difference 
approximation).  The  truncation  error  for  this  approximation  is 

T(tn)  =  At2  (^£(4)  +  +  l-E ^  . 

Therefore,  since  the  semi-discrete  form  is  0(h2),  this  approximation  method  overall 
is  0(h2)  when  At  =  O(h). 


3.2.2  Lorentz  Media 


For  Lorentz  media  we  convert  (3.5a)-(3.5b)  into  a  first  order  system  of  equations  in  three 
unknowns,  X  =  [e,p,  d,  q]T,  where  d  =  e,  q  =  p,  resulting  in 

MX  +  KX  =  J,  (3.19) 


with 

'  I  0  0  O' 

0/00 

M  = 

0  0  Coo  M  0 

0  0  0  / 


(3.20) 
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0 

0 

-i 

0 

0 

0 

0 

-i 

(c2K  +  edu)%M) 

,  ,2 

—-M 

£0 

,  (3.21) 

— 

B 

0 

A 

J=[  0  0 

0 

eo 

r 

(3.22) 

As  before,  we  apply  a  theta-scheme  to  MX  +  KX  =  J  to  obtain 

(M  +  9AtK)Xn+1  =  (M  -  (1  -  9)AtK)Xn  +  ( 9Jn+1  +  (1  -  6)  Jn),  (3.23) 

For  9  =  0.5  the  scheme  can  be  written  as 

(M  +  ^AT)An+1  =  (M  -  yR)Xh  +  l-{  Jn+l  +  Jn).  (3.24) 

where,  Jn  =  J(nAt)  and  A"  =  X(nAt). 

Again  we  are  assuming  a  fixed  At,  therefore  we  solve  for  the  LU  factorization  at  the 
beginning  and  use  back  substitution  at  each  time  step. 

4  Stability  Analysis 

Fourier  analysis  is  an  important  tool  in  the  study  of  stability  of  finite  difference  and  finite 
element  schemes.  The  Fourier  transform  of  a  function  gives  an  alternative  representation  of 
the  function,  and  one  can  infer  certain  properties  of  a  function  from  its  Fourier  transform. 
The  Fourier  inversion  formula  represents  a  function  as  a  superposition  of  waves  eluJZ  with 
different  amplitudes  that  are  given  by  the  Fourier  transform.  Under  the  Fourier  transform 
the  operation  of  differentiation  is  converted  into  the  operation  of  multiplication  by  ica. 

An  important  application  of  Fourier  analysis  is  the  von  Neumann  analysis  of  stability 
of  difference  schemes.  With  the  use  of  Fourier  analysis  we  can  give  necessary  and  sufficient 
conditions  for  the  stability  of  these  schemes.  For  a  difference  scheme,  advancing  the  solu¬ 
tion  of  the  scheme  by  one  time  step  is  equivalent  to  multiplying  the  Fourier  transform  of 
the  solution  by  an  amplification  factor.  The  amplification  factor  is  so  called  because  its 
magnitude  is  the  amount  that  the  amplitude  of  each  frequency  in  the  solution,  given  by  the 
Fourier  transform  of  the  solution  at  time  step  n,  is  amplified  in  advancing  the  solution  to 
time  step  n  +  1.  All  the  information  of  a  scheme  is  contained  in  its  amplification  factor. 
In  particular,  the  stability  and  accuracy  of  schemes  can  be  determined  from  the  amplifica¬ 
tion  factor.  Using  Fourier  analysis  on  the  scheme  to  calculate  the  amplification  factor  f  is 
equivalent  to  relacing  discrete  values  of  any  unknown  field  vector  v^,  at  time  t  =  nAt  and 
spatial  position  z  =  mAz,  in  the  scheme  by  (fneime  for  each  value  of  n  and  m.  The  resulting 
equation  is  a  polynomial  in  f  and  is  called  the  stability  polynomial  or  the  characteristic 
polynomial  for  the  scheme.  The  stability  polynomial  usually  depends  on  the  discretization 
parameters,  such  as  the  time  step  and  the  mesh  step  size  as  well  as  the  medium  parameters. 
The  roots  of  this  polynomial  can  be  obtained  and  will  determine  the  stability  of  the  scheme. 
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Since  the  roots  of  the  polynomial  are  the  amplification  factors  of  the  scheme,  the  scheme 
will  be  stable  if  the  amplitude  of  the  roots  is  less  than  or  equal  to  1. 

We  follow  [28]  to  derive  stability  results  for  the  finite  element  schemes  presented  here. 
Since  the  boundary  conditions  do  not  affect  the  stability  and  dispersion  properties  of  the 
scheme  in  the  interior  of  the  domain,  we  neglect  the  effects  of  boundary  conditions  (B  =  0 
in  (3.1))  in  our  analysis.  As  the  stability  and  dispersion  properties  of  the  scheme  are  also 
independent  of  the  source  Js,  we  take  J  =  0  in  (3.1)  without  loss  of  generality. 


4.1  Free  Space 

For  the  stability  analysis  in  free  space  we  have  P  =  0.  Thus,  in  equation  (3.1)  we  also 
substitute  P  =  0.  In  free  space  we  have  the  equation 

HoeoMe  +  Ke  =  0.  (4-1) 

In  the  above  we  have  also  assumed  that  =  1.  Using  c  =  1/^/ioMo,  and  rewriting  (4.1)  in 
first  order  form  in  the  variables  X  =  [e,  d]T1  where  d  =  e,  we  have 


/  0 

o  Am 


x  + 


0  -I 
I\  0 


X  =  0. 


Applying  a  theta  scheme  to  (4.2)  we  get 


I  -0A  tl 

QAtK  \M 

Cz 


Xn+1 


I  (1  —  9)  At  I 

-(1  -  6)AtK  ± M 


Xn  =  0. 


(4.2) 


(4.3) 


We  note  the  following  identities  associated  with  the  application  of  the  mass  and  the 
stiffness  matrices  on  vectors  ft™  =  ft(nel^k^Az\ 

1.  For  the  mass  matrix  M  we  have 


Mq 


Az 


>1-1 


—ikAz 


2A^ 

A^ 

+  3  <?  + 

- ( 

6 

+  4  +  eifcAz) 

Az 

6 

r 


i+i 


A, 

(3  —  2  sin2  (kAz/2))^4>r- 
3  J 


(4.4) 


2.  For  the  stiffness  matrix  K  we  have 


K‘i>j 


+ 4-j?  - 


n  jji 

A  z 
-1 


(e-ikAz  -2  +  eikAz^j  —ft 


A  z 


sin2  (kAz/2))ft 


(4.5) 
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We  also  define  two  quantities 

K  =  3-2sin2^^V  (4.6a) 

r]  =  3^2sin2  ,  (4.6b) 

where  the  Courant  number  is 

y  =  (e«  =  i).  (4.7) 

In  a  general  dispersive  material  v  =  In  air  =  1.  To  determine  the  stability 

criterion  for  the  finite  element  scheme  in  free  space  we  now  substitute 


X?  = 


1 

S-o 

1 _ 

i - 

dn 

L  aj  J 

d 

n  i(kjAz) 


Cne 


(4.8) 


into  the  discrete  equations  (4.3)  with  8  =  1/2.  Here  k  is  the  wave  number,  Using  the 
identities  (4.4)  and  (4.5)  we  obtain  a  homogeneous  linear  system  of  the  form 


AX  =  0, 


(4.9) 


with 


A  = 


2A  t 
~Kz 


sin 


C-1 

^)(C  +  D 


-f(C  +  i) 


(3-2*i“2(^))  (<-!)_ 


Az 

5?|3-2sm 


(4.10) 


and 


X  = 


e,  d 


-l  T 


By  setting  the  determinant  of  A  to  be  zero,  we  obtain  the  stability  polynomial 


(4.11) 


C  -2C 


4^'i+i=o. 

+  TJZ 


(4.12) 


From  (4.12)  we  can  show  that  |£|  =  1  always,  regardless  of  the  medium  parameters.  This 
implies  that  the  finite  element  scheme  with  the  theta  method  (8  =  1/2)  in  free  space  is 
unconditionally  stable  as  well  as  non-dissipative. 


4.2  Debye  Media 

To  determine  the  stability  conditions  for  the  finite  element  scheme  described  as  Method  1 
in  Section  3.2.1  for  Debye  media  we  substitute 


-  pn  ~ 

^3 

i - 

i _ 

II 

k*  "c~'5 

i 

_ i 

- 1 

_ i 

Qn^ikjAz 


(4.13) 
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in  the  discrete  equations  (3.13)  (with  J  =  0).  As  in  the  case  of  free  space,  we  obtain  a 
homogeneous  system  of  the  type  AX  =  0.  We  then  set  det(A)  =  0  to  obtain  the  stability 
polynomial 


03C3  +  CL2C 2  +  a\  C  +  ao  —  0,  (4.14) 

where  the  coefficients  of  the  stability  polynomial  are  given  as 

03  =  f]~{hT  +  2)  +  n2(hTes  +  2eoo),  (4.15) 

02  =  i]2(3hT  +  2)  -  K2(hTes  +  6eoo),  (4.16) 

01  =  r/2(3 hT  -  2)  -  K2(hTes  -  6eoo),  (4.17) 

a0  =  r]2(hT  -  2)  +  n2(hTes  -  2eoo),  (4.18) 


where  77,  and  k  are  as  defined  in  (4.6b),  and  (4.6a),  respectively,  and  hT  =  A t/r. 

To  determine  the  stability  polynomial  for  the  finite  element  scheme  described  in  Method 
2  in  Section  (3.2.1),  we  substitute  (4.13)  in  the  discrete  equations  (3.17).  Following  the 
procedure  discussed  above  we  obtain  the  stability  polynomial 


63C3  +  62C2  +  &iC  +  &0  =  0,  (4.19) 

with  coefficients 


(hT  +  2)  +  k  (hTes  +  2eoo)  ^  , 

(4.20) 

2  13 

r)2{3hT  +  2)  K~(hTes  +  6eoo)  +  "  ^  , 

(4.21) 

K2h3fj 

r]2(3hT  2)  K-(hTes  6eoo)  +  *  , 

(4.22) 

r, f(hT  2)  +  K2(hTes  26^)  *  , 

(4.23) 

If  we  neglect  terms  involving  h 3,  the  stability  polynomial  of  Method  2  for  Debye  is  the 
same  as  in  Method  1.  Thus  for  small  values  of  hT,  the  two  methods  have  the  same  stability 
properties.  In  the  rest  of  this  section  we  refer  to  Method  1  as  the  finite  element  scheme 
and  we  will  compare  stability  properties  of  this  scheme  with  those  of  the  finite  difference 
scheme  [18]  analyzed  in  [28]. 

In  Figure  1  we  plot  the  absolute  value  of  the  largest  root  of  (4.14)  for  the  finite  element 
scheme  with  v  =  1  as  a  function  of  kAz.  In  Figure  2  we  plot  the  the  absolute  value  of  the 
largest  root  of  the  stability  polynomial  of  the  finite  difference  scheme  JHT  [18],  with  v  =  1 
as  a  function  of  kAz. 

In  Figure  3  (left)  we  plot  the  absolute  value  of  the  largest  root  of  (4.14)  for  the  finite 
element  scheme  with  v  =  1  as  a  function  of  the  wave  number  k.  In  Figure  3  (right)  we  plot 
the  the  absolute  value  of  the  largest  root  of  the  stability  polynomial  of  the  finite  difference 
scheme  JHT  [18],  with  v  =  1  as  a  function  of  k. 

To  generate  these  plots  we  assumed  the  following  values  of  the  physical  parameters: 
r  =  8.1  x  10~12,  Coo  =  1,  and  es  =  78.2.  The  time  step  At  is  determined  by  the  choice  of 
hT  and  the  physical  parameter  r.  These  plots  show  varying  values  of  hT  from  0.1  to  0.001. 
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Debye  stability  for  FEM  with  v=1  Debye  stability  for  FEM  with  v=1 


Figure  1:  Plot  of  max|£|  versus  kAz  for  hT  =  {0.1,  0.01,  0.001}  for  the  finite  element  scheme 
in  a  Debye  medium,  with  v  =  1.  (Right  plot  is  zoom  of  left  plot.) 


From  the  plots  we  can  see  that  the  dissipation  of  the  numerical  schemes  can  be  reduced 
by  decreasing  hT.  For  stability  and  least  dissipation,  hT  =  0.001  is  recommended.  For 
the  finite  element  scheme  increasing  u  from  1  to  16  does  not  seem  to  change  the  stability 
behaviour  of  the  scheme.  This  suggests  the  unconditional  stability  of  the  finite  element 
scheme.  However,  the  finite  difference  scheme  is  conditionally  stable  and  has  the  stability 
criteria  u  <  1.  The  stability  criteria  for  the  finite  difference  scheme  have  been  derived  in 
[28,  10], 

In  Figure  4  we  plot  max|£|  versus  kAz  for  hT  =  0.1  (left)  and  hT  =  0.3  for  the  finite 
element  schemes  of  Method  1  (FEM1)  and  Method  2  (FEM2)  in  a  Debye  medium,  with 
v  =  1.  From  this  plot  we  can  see  a  slight  difference  in  the  two  schemes  when  hT  =  0.3. 
However,  for  hT  =  0.1  the  two  schemes  produce  the  same  plots. 
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Debye  stability  for  FD  with  v=1 


Debye  stability  for  FD  with  v=1 


Figure  2:  Plot  of  max|£|  versus  kA.z  for  hT  =  {0.1,0.01,0.001}  for  the  finite  difference 
scheme  JHT  in  a  Debye  medium,  with  v  =  1.  (Right  plot  is  zoom  of  left  plot.) 


Debye  stability  for  FEM  with  v=1  Debye  stability  for  FD  with  v=1 


Figure  3:  Plot  of  max|<(|  versus  the  wave  number  k  for  hT  =  {0.1,0.01,0.001}  for  the  finite 
element  scheme  (left)  and  the  finite  difference  scheme  (right)  in  a  Debye  medium,  with 
v  =  1. 
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Debye  stability  comparison  with  v=1 ,  h=0.3  Debye  stability  comparison  with  v=1 ,  h=0.1 


Figure  4:  Plots  of  max|£|  versus  kAz  for  hT  =  0.1  (left)  and  hT  =  0.3  for  the  finite  element 
schemes  of  Method  1  (FEM1)  and  Method  2  (FEM2)  in  a  Debye  medium,  with  u  =  l. 


4.3  Lorentz  Media 

To  determine  the  stability  conditions  for  Lorentz  media  we  substitute 


XJ  = 


1 

<13 

1 _ 

e 

P? 

P 

dU; 

(1 

J 

l^\ 

.  Q  . 

nQikjAz 


(4.24) 


in  the  discrete  equations  (3.24).  As  in  the  case  of  free  space,  we  obtain  a  homogeneous 
system  of  the  type  AX  =  0.  We  then  set  det(*4)  =  0  to  obtain  the  stability  polynomial 


&4:C^  T  o-sC^  T  ®2C2  T  aiC  +  no  —  0,  (4.25) 

where  the  coefficients  of  the  stability  polynomial  are  given  as 

04  =  tT(27t2/iq  +  hT  +  2)  +  K~{2-K~hles  +  hTe  oo  +  2600)  (4.26) 

a3  =  r]2{8ir2hl  +  2 hT)  -  K2eoo( 8  +  2 hT)  (4.27) 

a2  =  r/2(127r2/rQ  —  4)  —  K2(47r2/iQes  —  12eoo)  (4.28) 

CLi  =  r?2(87 T2h20  -  2 hT)  -  K2eoo(8  -  2 hT)  (4.29) 

a0  =  r/2(27r2/ig  -  hT  +  2)  -  K2(2ir2h%es  -  hre^  +  26^),  (4.30) 


where  77,  and  k  are  as  defined  in  (4.6b),  and  (4.6a),  respectively,  hT  =  A t/r,  and  ho  =  At/To, 
where  To  =  2n/ioo- 

We  plot  the  absolute  value  of  the  largest  root  of  (4.25)  for  v  =  1  versus  kAz  for  the 
finite  element  scheme  in  Figures  5,  and  for  the  finite  difference  schemes  JHT  [18]  and  KF 
[19]  in  Figures  6  and  7. 
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Lorentz  stability  for  FEM  with  v=1  Lorentz  stability  for  FEM  with  v=1 


Figure  5:  Plot  of  max|£|  versus  kAz  for  hr  =  {0.1,  0.01,  0.001}  for  the  finite  element  scheme 
in  a  Lorentz  medium,  with  v  =  1.  (Right  plot  is  zoom  of  left  plot.) 


We  then  plot  the  absolute  value  of  the  largest  root  of  (4.25)  for  v  =  1  versus  the  wave 
number  k  for  the  finite  element  scheme  in  Figures  8,  and  for  the  finite  difference  schemes 
JHT  and  KF  in  Figures  9. 

To  generate  these  plots  we  assumed  the  following  values  for  physical  parameters:  r  = 
1.786  x  10~16,  Coo  =  1,  es  =  2.25,  and  uiq  =  4  x  1016.  For  the  Lorentz  medium,  all  time 
scales  must  be  properly  resolved,  therefore  the  time  step  At  is  determined  by  the  choice  of 
either  hT  or  ho,  whichever  is  most  restrictive.  For  the  current  parameter  values,  Tq  <  t, 
thus  ho  is  used.  The  plots  show  varying  values  of  ho  from  0.1  to  0.001.  From  the  plots  we 
can  see  that  the  dissipation  of  the  numerical  schemes  can  be  reduced  by  decreasing  ho-  For 
stability  and  least  dissipation,  ho  =  0.01  is  recommended. 

As  in  the  case  of  the  finite  element  method  for  Debye  media,  we  see  that  increasing  v 
from  1  to  16  does  not  affect  the  stability  properties  of  the  finite  element  scheme  for  Lorentz 
media.  This  suggests  the  unconditional  stability  of  the  finite  element  method.  However  the 
finite  difference  scheme  is  again  conditionally  stable  with  the  criteria  v  <  1.  The  stability 
criteria  for  the  finite  difference  scheme  for  Lorentz  media  have  also  been  derived  in  [28,  10]. 
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Lorentz  stability  for  JHT  with  v=1  Lorentz  stability  for  JHT  with  v=1 


Figure  6:  Plot  of  max|£|  versus  kAz  for  hT  =  {0.1,  0.01,  0.001}  for  the  JHT  finite  difference 
scheme  in  a  Lorentz  medium,  with  v  =  1.  (Right  plot  is  zoom  of  left  plot.) 


Lorentz  stability  for  KF  with  v=1  Lorentz  stability  for  KF  with  v=1 


Figure  7:  Plot  of  max|£|  versus  kAz  for  hT  =  {0.1,0.01,0.001}  for  the  KF  finite  difference 
scheme  in  a  Lorentz  medium,  with  v  =  1.  (Right  plot  is  zoom  of  left  plot.) 
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Lorentz  stability  for  FEM  with  v=1 


Figure  8:  Plot  of  max|£|  versus  the  wave  number  k  for  hT  =  {0.1,0.01,0.001}  for  the  KF 
finite  element  scheme  in  a  Lorentz  medium,  with  v  =  1. 


Lorentz  stability  for  JHT  with  v=1  Lorentz  stability  for  KF  with  v=1 


Figure  9:  Plot  of  max|£|  versus  the  wave  number  k  for  hT  =  {0.1,0.01,0.001}  for  the  JHT 
(left)  and  KF  (right)  finite  difference  schemes  for  Lorentz  medium,  with  v  =  1. 
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5  Analysis  of  Dispersion  and  Phase  Error 


The  numerical  approximation  of  time-dependent  wave  problems  introduces  errors  which 
involve  dissipation,  dispersion,  and  anisotropy.  The  attenuation  of  the  amplitude  of  the 
plane  wave  is  referred  to  as  dissipation.  The  numerical  model  produces  waves  that  propa¬ 
gate  at  incorrect  speeds.  The  dependence  of  the  velocity  of  propagation  of  the  numerical 
sinusoidal  waves  on  frequency  is  termed  as  dispersion,  while  the  dependence  of  the  velocity 
upon  direction  is  referred  to  as  numerical  anisotropy. 

A  dispersive  model  admits  plane  wave  solutions  of  the  form  for  wj1jcj1  the 

speed  of  propagation,  governed  by  the  wave  number  k,  is  not  independent  of  the  frequency 
u).  For  time-harmonic  waves,  numerical  dispersion  results  in  the  creation  of  a  numerical 
phase  velocity  error,  or  phase  error,  in  the  solution.  This  is  due  to  the  incorrect  modeling 
of  the  sinusoidal  behavior  of  the  propagating  wave,  for  example,  the  piecewise  polynomial 
approximation  of  a  finite  element  method  does  not  exactly  match  a  sine  or  cosine  function. 
Dispersion  is  present  in  numerical  approximation  methods  such  as  finite  difference/finite 
element  methods  even  in  the  absence  of  any  dispersion  in  the  actual  media.  Errors  are 
cumulative,  thus  as  waves  propagate  over  long  distances  the  numerical  solution  becomes 
corrupted  and  may  completely  deviate  from  the  correct  solution. 


5.0.1  Free  Space 

Substituting  a  solution  of  the  form 

e(t,  z )  =  e^kz~^ 

into  equation  (2.8)  results  in  the  dispersion  relation  of  the  continuous  model  in  free  space 
being 

^ex(w)  =  (5-1) 


where  denotes  the  wavenumber  in  free  space  (where  V  in  the  superscript  denotes 
“vacuum” )  for  the  exact  equations  (EX  in  the  subscript  denotes  “exact” ) .  To  determine  the 
dispersion  relation  for  the  discretized  model  using  the  finite  element  method,  we  substitute 


\rn _ 

^3  ~ 


\  e"l 

e 

L  aj  J 

d 

Az—umAt) 


(5.2) 


where  k/\  is  the  numerical  wave  number,  into  the  discrete  equation  (4.3).  Using  the  identities 
(4.4),  (4.5)  as  well  as  the  following  two  trignometric  identities 

g— iajAt/2  +  ekjAi/2  =  2cos(wAt/2),  (5.3a) 

e— icAt/2  _  ei^Ai/2  =  _2isin(cuAf/2).  (5.3b) 

we  obtain  the  linear  system 

AX  =  0,  (5.4) 

with 

2isin(^)  At  cos  (^) 

A  A/\+  _  /  Jc a  A?\  /tuAi\  A z 


4 At  2  ( k/\Az\ 

— —  sin  -  cos 

Az  V  2  / 


V  2 


-^!3-2sin 


.^2  ( &aA z\ 


\  2 


2isin 


u;At\ 
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and 


X=[e,d  ]' 


(5.5) 


Setting  the  determinant  of  A  =  0  we  obtain  a  relation  between  kA  and  u,  which  is  the 
numerical  dispersion  relation  for  the  finite  element  scheme  in  free  space.  This  relation  is 


2  ,  .  .  .  / 2  u2  cos2  (wAi/2)\ 

v  1  V3  sm2  (coAt/2)  ) 


-i 


Solving  for  kA  in  the  above  we  get 


(5.6) 


kA  =  kp e(w)  =  —  sin  1 


1 


1 2  v2  cos2  (ioAt/2) 


\  y  3  sin2  (u At/ 2)  ) 

The  dispersion  relation  for  the  FDTD  scheme  is  given  to  be 

sin  (uAt/2) 


which  implies  that 


sin  (. kAAz/2 )  = 


kA  =  kpD(uj)  =  —  sin  1 


sin  (ujAt/2) 


We  will  define  the  phase  error  as 


<f>(u;Af)  = 


&Ex(wAf)  —  kA(ujAt) 


kpx^At) 


(5.7) 


(5.8) 

(5.9) 

(5.10) 


where  for  the  finite  element  scheme  in  a  vacuum,  we  have  kA  =  kpE  as  defined  in  (5.7), 
whereas  for  the  FDTD  schemes  we  have  kA  =  kpD  as  defined  in  (5.9). 

In  Figures  10  and  11  we  plot  the  phase  error  in  a  vacuum  for  the  finite  element 
scheme  with  the  theta  method  {v  =  y^l/2, 1,4, 16)  and  for  the  FDTD  scheme  {y  = 
0.6,  1/2,  y^/3, 1),  respectively. 

To  see  why  v  ~  .7  has  the  least  dispersion  for  FEM  (see  plot  on  the  right  in  Figure  10) 
and  why  v  ~  1  has  the  least  dispersion  for  FDTD  (see  plot  on  the  right  in  Figure  11),  it  is 
helpful  to  plot  the  relations  in  equations  (5.6)  and  (5.8)  versus  the  continuous  model  values. 
We  define 


(5.11) 


where  we  have  substituted  kp-^Az  =  vujAt.  We  similarly  define  7pE  =  sin2  (kppAz/2) 
and  7p£,  =  sin2(fcEDAz/2)  using  the  definitions  in  (5.7)  and  (5.9),  respectively.  Figure  12 
displays  plots  each  of  these  72(cjA t)  functions  for  various  values  of  v.  For  the  continuous 
model  (left  plot),  v  has  the  effect  of  moving  the  location  of  the  maximum  value,  y2  =  1. 
For  the  finite  difference  case  (right  plot)  the  location  of  the  maximum  does  not  change, 
although  the  value  of  the  maximum  does.  For  v  =  1  the  curve  coincides  exactly  with  the 
continuous  case.  For  the  finite  element  method  (middle  plot)  the  location  of  the  maximum 
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to  At  to  At 


Figure  10:  Plot  of  the  phase  error  <t>  versus  u>At  for  the  finite  element  scheme  in  freespace 
with  v  =  {-y/l/2, 1,4, 16}.  For  v  =  yH/2,  FEM  has  the  least  dispersion.  (Right  plot  is 
zoom  of  left  plot.) 


to  At  to  At 


Figure  11:  Plot  of  the  phase  error  <f>  versus  uvAt  for  the  finite  difference  scheme  in  freespace 
with  v  =  {1,  -y/ 2/3,  \J  1/2,  .6}.  For  v  =  1  FDTD  has  zero  dispersion.  (Right  plot  is  zoom 
of  left  plot.) 


does  not  change,  nor  does  the  value.  However,  as  this  value  is  fixed  at  1.5,  it  will  never 
coincide  exactly  with  the  continuous  case  for  any  value  of  v. 

To  determine  the  best  value  of  v  for  the  finite  element  method,  we  first  note  that  7^x 
can  be  expanded  as 


1/4 


U)‘ 


1 At 2 


-  1/48 


U! 


AAt 4 


+ 


U! 


; At 6 


1440 


+  0  (u)8At8)  . 


(5.12) 


As  both  the  expansions  of  7fd  and  7fe  match  up  to  the  second  order  coefficient,  it  is  the 
fourth  order  Taylor  coefficient  that  determines  how  well  the  discretization  method  matches 
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Figure  12:  Plot  of  the  72  versus  uiAt  for  the  continuous  model,  the  finite  elment  scheme, 
and  the  finite  difference  scheme  (left  to  right). 


the  continuous  model.  For  the  finite  difference  method  we  have 

cED  =  —1/48  (5.13) 

whereas  for  the  finite  element  method  we  have 

<T  =  l/24=^.  (5.14) 

For  v  =  1,  D  =  cfx  =  1/48,  but  cEE  =  0.  If  we  solve  (cfx  —  cfE)(^)  =  0  we  find 
that  v  =  y/l/2  is  the  value  for  which  the  finite  element  method  most  closely  matches  the 
continuous  model. 

In  Figure  13  we  plot  the  fourth-order  Taylor  coefficient  of  sin(fcAz/2)  as  a  function  of  v 
for  k  =  {A:EX,  &pE,  A:pD},  i.e. ,  the  k  values  for  the  continuous  model,  the  finite  element  scheme 
and  finite  difference  scheme,  respectively.  As  mentioned  above,  for  v  =  y^l/2  the  coefficient 
corresponding  to  the  finite  element  scheme  is  equal  to  that  of  the  value  corresponding  to  the 
exact  k.  This  suggests  that  least  dispersion  will  be  achieved  in  the  finite  element  scheme 
when  v  =  ^/l/2.  For  v  =  1  the  coefficient  corresponding  to  the  finite  difference  scheme 
is  equal  to  that  of  the  value  corresponding  to  &EX.  It  is  known  that  this  is  the  value  of 
v  for  which  finite  difference  schemes  exhibit  least  dispersion.  Further,  for  v  =  -y/2/3  the 
coefficient  corresponding  to  the  finite  element  scheme  is  equal  to  that  of  the  finite  difference 
scheme.  The  dispersion  curves  are  nearly  identical  for  the  two  schemes  when  v  =  \J 2/?>. 
The  plots  of  q2  for  these  values  of  v  are  displayed  in  Figure  14. 

In  Figures  15  and  16  we  plot  the  phase  velocity  vp  (scaled  by  p)  versus  u>A t  for  the 
finite  element  scheme  with  the  theta  method  and  the  phase  error  for  the  FDTD  schemes, 
respectively. 
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c4 


Figure  13:  Value  of  the  4th  order  Taylor  coefficient  of  sin2(^^)  for  k  =  {fcgX,  ^fe>  &fD}- 


Figure  14:  Plot  of  the  y2  versus  uiAt  for  the  continuous  model,  the  finite  elrnent  scheme, 
and  the  finite  difference  scheme  (Left  plot  uses  v  =  1 ,  center  plot  uses  v  =  \J 2/3,  right  plot 
uses  v  =  V/I72)- 
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Free  Space  velocity  for  FEM  Free  Space  velocity  for  FEM 


Figure  15:  Plot  of  the  velocity  vp  (scaled  by  -)  versus  uAt  for  the  finite  element  scheme  in 
freespace  with  v  =  {1/1/2, 1,4, 16}.  For  v  =  1/2,  FEM  has  the  velocity  closest  to  one. 

Note  that  for  v  <  i/l/2,  the  velocity  vv  is  greater  than  c. 


Free  Space  velocity  for  FD  Free  Space  velocity  for  FD 


Figure  16:  Plot  of  the  velocity  vp  (scaled  by  versus  wA t  for  the  finite  difference  scheme 
in  freespace  with  v  =  {1,  yj 2/3 ,  y/l/2,  .6}.  For  v  =  1  FDTD  always  has  velocity  one. 
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5.0.2  Debye  Media 

The  dispersion  relation  for  the  continuous  Debye  model  is  given  by 


(5.15) 


where 

M  =  6^A  •  ^ ,  (5-16) 

A  — 

is  the  relative  complex  permittivity  of  the  Debye  medium,  A  =  1/r.  To  determine  the 
numerical  dispersion  relation  for  the  finite  element  method  described  as  Method  1  in  Section 
3.2.1  applied  to  a  Debye  media,  we  substitute 


\eU 

e 

x?  = 

Pj 

= 

P 

U  J 

.  d  _ 

Az—umAt) 


(5.17) 


in  the  discrete  equations  (3.13)  (with  J  =  0).  We  assume  again  that  Coo  =  1.  As  in  the 
case  of  free  space,  we  obtain  a  homogeneous  system  of  the  type  AX  =  0.  We  then  set 
det(*4)  =  0  to  obtain  a  relation  between  the  numerical  wavenumber  and  the  frequency 
ui.  Solving  for  &a  hr  this  relation,  and  writing  in  a  form  comparable  to  (5.15),  we  obtain 
the  numerical  dispersion  relation  for  the  FEM1  scheme  applied  to  Debye  media  to  be 


kA  =  kpE(cj)  =  Az  Sin  1 


uja  Az 
~cT^~ 


-D 

'r,FE 


with  the  discrete  relative  complex  permittivity  given  by 

£s,  A^A  ~  itUA 


,D  _ 

er,FE  — 


Aa  —  iu>A 


where  the  discrete  medium  parameters  are 


and 


(5.18) 


(5.19) 


es 

£sA  ~  a2/32 

(5.20) 

Aa  =  \cos(cuAt/2)/32a3 

(5.21) 

CO  A  =  WSujQ 

(5.22) 

sin(o;Ai/2) 

(5.23) 

ujAt/2 

a  = 

/2sin2(a;At/2)  2  A  lA 

f - }  +  cos2(caAt/2) J 

(5.24) 

P  = 

( 2es  sin2(cjAf/2)  2  \  1//2 

f  3^2  '  ’  +  cos2  {uAt/2)\  . 

(5.25) 

To  determine  the  numerical  dispersion  relation  for  the  finite  element  scheme  described  in 
Method  2  in  Section  (3.2.1),  we  substitute  (4.13)  in  the  discrete  equations  (3.17).  Following 
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the  procedure  discussed  above  we  obtain  the  numerical  dispersion  relation  for  the  FEM1 
scheme  applied  to  Debye  media  to  be 


sin2(/cAAz/2)  =  a 2 


2sin(?y)i  —  Acos(r/)esAt  +  cos(7/)e<i 
2sin(?y)i  —  Xcos(r])esAt/32a2  +  bh^  cos(r])eda2 


(5.26) 


where  r/  =  =  es  —  1  and  /iT  =  At/r.  If  we  neglect  terms  in  h^.,  then  the  expression 

(5.26)  reduces  to  (5.18).  Thus,  for  small  hT,  both  the  finite  element  methods  have  the  same 
numerical  dispersion  relations.  From  the  section  on  stability  analysis  we  have  seen  that  for 
low  dissapation  ht/m  needs  to  be  about  0.001  for  Debye  media  and  0.01  for  Lorentz  media. 
For  these  values  of  hT  both  finite  element  schemes  produce  the  same  dispersion  graphs. 

Both  the  FEM  schemes  misrepresent  the  continuous  model  parameters  A  and  es  dis¬ 
cretely  as  Aa  and  cS)a,  and  misrepresent  the  frequency  u  as  wa-  We  note  that  as  v 
increases  [y  >  y/ej),the  product  a(3  — ►  1,  and  eS)A  — ►  es-  The  discrete  parameter  Aa 
for  the  finite  element  method  is  a  function  of  the  continuous  model  parameter  es  via  the 
quantity  (3.  However,  for  the  regime  of  interest,  namely  uAt  small,  and  when  v  is  large 
(u  >  y TJ),  2es sin2(cuAf/2)/3^2  is  dominated  by  cos2 (co At/ 2).  Also  the  product  a/3  — ►  1, 
a  — >  1/ cos(cuA/2)  and  thus  Aa  — >  A.  Finally  for  u>At  small  and  v  large,  s^a  «  1  and 
thus  cua  — >  oj.  Hence,  the  choice  of  the  Courant  number  v  is  important  in  maintaining  low 
dispersion  error  in  the  discrete  FEM. 

We  compare  the  FEM  scheme  for  Debye  media  with  the  finite  difference  scheme  pre¬ 
sented  in  [18]  and  analyzed  in  [28]  (for  the  Debye  media,  this  method  is  again  equivalent 
to  the  scheme  in  [20]).  The  numerical  dispersion  relation  for  this  finite  difference  scheme 
applied  to  Debye  media  is  given  to  be 


A:  a  —  A/p£>(o;) 


Az 


sin 


-l 


lo  Az  nz 

-^VerTD  . 


(5.27) 


er,FD 


€sXa  ~  A 
Aa  —  iuj  a 


(5.28) 


If  written  in  the  form  of  (5.18),  this  would  correspond  to  the  following  discrete  representa¬ 
tions  of  the  continuous  model  parameters: 


es,  A  =  es  (5.29) 

Aa  =  Acos(cuAf/2),  (5.30) 

and  a  representation  of  the  frequency  by 

u;a  =  (5.31) 

We  compare  the  phase  error  for  the  finite  element  scheme  applied  to  Debye  media  to  phase 
error  for  the  finite  difference  scheme.  The  phase  error  is  plotted  against  values  of  cuAt  in 
the  range  [0,7r].  We  note  that  uiAt  =  27r/AipPp,  where  AVPPP  is  the  number  of  points  per 
period,  and  is  related  to  the  number  of  points  per  wavelength  lVppw  via 

ArPPw  =  ^  A  pPP .  (5.32) 
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Thus,  for  v  <  1,  the  number  of  points  per  wavelength  is  always  less  than  or  equal  to  the 
number  of  points  per  period,  and  conversely  for  v  >  1.  Note  that  the  number  of  points  per 
wavelength  in  the  range  [7t/4,7t]  is  8  to  2  points  per  period.  We  are  more  interested  in  the 
range  [0,  vr/4]  which  involves  more  than  8  points  per  period  (or  equvivalently  more  than  8 
points  per  wavelength). 

To  generate  the  plots  below  we  have  used  the  following  values  for  the  medium  parame¬ 
ters: 


eoo  =  1  (5.33) 

es  =  78.2  (5.34) 

r  =  8.1  x  10”12  sec.  (5.35) 

Figures  17-19  plot  the  phase  error  $  versus  uAt  for  the  FEM  scheme  applied  to  the 
Debye  model,  with  various  values  of  hT  and  u.  Figures  20-22  plot  the  same  for  the  finite 
difference  method. 

Figures  23-25  plot  the  phase  error  4>  versus  the  frequency  u>  for  the  FEM  scheme  applied 
to  the  Debye  model,  with  various  values  of  hT  and  v.  Figures  26-28  plot  the  same  for  the 
finite  difference  method.  In  these  plots  we  can  see  that  in  the  finite  element  scheme  the 
phase  error  reduces  as  v  increases,  even  beyond  1. 

In  Figure  31  plots  of  the  phase  error  in  the  two  finite  element  schemes  described  in 
Method  1  and  Method  2  in  Section  3.2.1  are  shown  with  v  =  1  and  hT  =  0.3  (left)  and 
hT  =  0.1  (right).  As  can  be  seen  in  these  plots  the  phase  error  in  both  the  schemes  is 
identical.  A  slight  difference  is  noticeable  in  the  case  hT  =  0.3,  however  for  hT  =  0.1  the 
two  graphs  are  identical. 

In  Figures  32-37  we  plot  the  real  and  imaginary  parts  of  the  relative  complex  permittiv¬ 
ity  for  Debye  media  for  the  continuous  equations  (exact  values),  the  finite  element  scheme 
with  v  =  {y/l/2, 1,4},  hT  =  {0.1,0.01},  and  finally  the  finite  difference  scheme  with  v  =  1 
fixed  and  hT  =  {0.1,0.01}.  For  hT  =  0.1,  as  v  is  increased  the  discrete  permittivites  of 
the  finite  element  scheme  approach  the  exact  values.  For  hT  =  0.01,  the  agreement  of  the 
discrete  permittivities  with  the  exact  values  is  better  than  with  hT  =  0.1,  for  each  value  of 
v.  For  v  =  4  we  see  the  best  agreement  of  the  discrete  real  and  imaginary  permittivities 
with  the  exact  values.  This  is  again  seen  in  the  plots  of  the  discrete  values  of  the  parameters 
A  and  es  in  Figures  39  and  38,  respectively.  In  both  these  figures  we  see  that  the  discrete 
parameters  in  the  finite  element  scheme  have  better  agreement  with  the  exact  values  as 
v  is  increased  and  hT  is  decreases.  However  the  discrete  value  of  the  frequency  tu,  Figure 
40,  seems  to  have  a  better  agreement  with  the  exact  value  when  v  =  1.  For  the  finite 
difference  scheme,  the  discrete  parameters  do  not  depend  on  the  value  of  zx  However  as  hT 
is  decreased  they  agree  better  with  the  exact  values.  For  the  particular  values  tested  here, 
it  seems  that  the  value  of  v  that  correctly  represents  A  will  sufficiently  model  the  complex 
permittivity. 
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Debye  dispersion  for  FEM  with  h  =0.1  Debye  dispersion  for  FEM  with  h  =0.1 


Figure  17:  Plot  of  the  phase  error  </>  versus  ujAt  for  the  finite  element  scheme  of  the  Debye 
model  with  v  =  {yJl/2, 1,4, 16}  using  hT  =  0.1.  (Right  plot  is  zoom  of  left  plot.) 


Debye  dispersion  for  FEM  with  h  =0.01  Debye  dispersion  for  FEM  with  h  =0.01 


Figure  18:  Plot  of  the  phase  error  </>  versus  uiAt  for  the  finite  element  scheme  of  the  Debye 
model  with  v  =  {y/l/2, 1,4, 16}  using  hT  =  0.01.  (Right  plot  is  zoom  of  left  plot.) 
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to  At  to  At 


Figure  19:  Plot  of  the  phase  error  </>  versus  ujAt  for  the  finite  element  scheme  of  the  Debye 
model  with  v  =  {yR/2, 1,4, 16}  using  hT  =  0.001.  (Right  plot  is  zoom  of  left  plot.) 


Debye  dispersion  for  FD  with  h  =0.1 


Debye  dispersion  for  FD  with  h  =0.1 


Figure  20:  Plot  of  the  phase  error  <j>  versus  ujAt  for  the  finite  difference  scheme  of  the  Debye 
model  with  v  =  {1,  -y/ 2/3,  -^/l/2,  .6}  using  hT  =  0.1.  (Right  plot  is  zoom  of  left  plot.) 
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Figure  21:  Plot  of  the  phase  error  </>  versus  ujAt  for  the  finite  difference  scheme  of  the  Debye 
model  with  v  =  {1,  -y/ 2/3,  -^/l/2,  .6}  using  hT  =  0.01.  (Right  plot  is  zoom  of  left  plot.) 


to  A  t  to  A  t 


Figure  22:  Plot  of  the  phase  error  </>  versus  uAt  for  the  finite  difference  scheme  of  the  Debye 
model  with  v  =  {1,  -y/ 2/3,  -^/l/2,  .6}  using  /iT  =  0.001.  (Right  plot  is  zoom  of  left  plot.) 
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Figure  23:  Plot  of  the  phase  error  cj)  versus  oj  for  the  finite  element  scheme  of  the  Debye 
model  with  hT  =  0.1  using  u  =  {-y/l/2, 1,4, 16}.  (Right  plot  is  log  of  left  plot.) 


Figure  24:  Plot  of  the  phase  error  cj)  versus  oj  for  the  finite  element  scheme  of  the  Debye 
model  with  hT  =  0.01  using  u  =  j-^/l/2, 1,4, 16}.  (Right  plot  is  log  of  left  plot.) 
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Figure  25:  Plot  of  the  phase  error  cj)  versus  oj  for  the  finite  element  scheme  of  the  Debye 
model  with  hT  =  0.001  using  v  =  {-y/l/2, 1,4, 16}.  (Right  plot  is  log  of  left  plot.) 


Debye  dispersion  for  FD  with  h  =0.1 


Debye  dispersion  for  FD  with  h  =0.1 


Figure  26:  Plot  of  the  phase  error  </>  versus  u>  for  the  finite  difference  scheme  of  the  Debye 
model  with  hT  =  0.1  using  v  =  {0.6,  -y/l/2,  ^2/3, 1}.  (Right  plot  is  log  of  left  plot.) 
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Figure  27:  Plot  of  the  phase  error  </>  versus  u>  for  the  finite  difference  scheme  of  the  Debye 
model  with  hT  =  0.01  using  u  =  {0.6,  \J  1/2,  -y/2/3, 1}.  (Right  plot  is  log  of  left  plot.) 


Figure  28:  Plot  of  the  phase  error  cj>  versus  u>  for  the  finite  difference  scheme  of  the  Debye 
model  with  hT  =  0.001  using  v  =  {0.6,  y^l/2,  ^2/3, 1}.  (Right  plot  is  log  of  left  plot.) 
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Figure  29:  Plot  of  the  phase  error  cj)  versus  uj  for  the  finite  element  scheme  of  the  Debye 
model  with  v  =  1  using  hT  =  {0.1,0.01,0.001}.  (Right  plot  is  log  of  left  plot.) 


Debye  dispersion  for  FD  with  v=1 


Figure  30:  Plot  of  the  phase  error  <jj  versus  uj  for  the  finite  difference  scheme  of  the  Debye 
model  with  v  =  1  using  hT  =  {0.1,0.01,0.001}.  (Right  plot  is  log  of  left  plot.) 
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Debye  dispersion  with  v=1 ,  h=0.3  Debye  dispersion  with  v=1 ,  h=0.1 


Figure  31:  Comparison  of  the  phase  error  in  the  two  finite  element  schemes  with  v  =  1  and 
hT  =  0.3  (left)  and  hT  =  0.1  (right). 
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Debye  with  =  0.1 ,  v  =  V(1/2) 


Figure  32:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Debye  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  y^l/2,  hT  = 
0.1  and  the  finite  difference  scheme  with  v  =  1  ,hT  =  0.1. 


Debye  with  h^  =  0.01 ,  v  =  V(1/2) 


Figure  33:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Debye  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  yD/2,  hT  = 
0.01  and  the  finite  difference  scheme  with  v  =  1  ,hT  =  0.01. 
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Debye  with  hx  =  0.1 ,  v  =  1 


Figure  34:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Debye  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  l,hT  =  0.1 
and  the  finite  difference  scheme  with  v  =  1,  hT  =  0.1. 


Debye  with  h^  =  0.01 ,  v  =  1 


Figure  35:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Debye  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  1,  hT  =  0.01 
and  the  finite  difference  scheme  with  v  =  l,hT  =  0.01. 
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Debye  with  hx  =  0.1 ,  v  =  4 


Figure  36:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Debye  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  4,  hT  =  0.1 
and  the  finite  difference  scheme  with  v  =  1,  hT  =  0.1. 


Debye  with  h^  =  0.01 ,  v  =  4 


Figure  37:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Debye  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  4,  hT  =  0.01 
and  the  finite  difference  scheme  with  v  =  1  ,hT  =  0.01. 
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Debye  with  h  =  0.1  Debye  with  h  =  0.01 
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Figure  38:  Plots  of  es  for  Debye  media  for  the  continuous  equations,  the  finite  element 
scheme  and  the  finite  difference  scheme  (u  =  1)  with  hT  =  0.1  (left)  and  hT  =  0.01  (right). 


Debye  with  h  =  0.1 


16 - FEM  v=V(1/2) 

- FEM  v=1 

FEM  v=4 


5  6  7 


Figure  39:  Plots  of  A  for  Debye  media  for  the  continuous  equations,  the  finite  element 
scheme  and  the  finite  difference  scheme  (u  =  1)  with  hT  =  0.1  (left)  and  hT  =  0.01  (right). 
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Figure  40:  Plots  of  the  frequency  u>  for  Debye  media  for  the  continuous  equations,  the  finite 
element  scheme  and  the  finite  difference  scheme  {y  =  1)  with  hT  =  0.1  (left)  and  hT  =  0.01 
(right). 


5.0.3  Lorentz  Media 


The  dispersion  relation  for  the  continuous  Lorentz  model  is  given  by 

^Ex(^)  =  —\J er  (w)' 

where  the  relative  complex  permittivity  for  Lorentz  media  is  given  to  be 

lo 2  —  esuiQ  +  iXu) 


er  (w)  = 


co2  —  ujn  +  iXuo  ’ 


(5.36) 


(5.37) 


To  determine  the  numerical  dispersion  relation  for  the  finite  element  method  applied  to  a 
Lorentz  medium  we  substitute 


yn  _ 
^3  ~ 


1 

1 _ 

e 

Pj 

P 

dr- 1 

d 

.  Q  . 

Az—cjnAt) 


(5.38) 


into  the  discrete  equations  (3. 19), (3. 20)  and  (3.21).  Again,  we  asssume  that  Coo  =  1  in  the 
following  analysis. 

As  before  we  obtain  a  homogeneous  system  of  the  type  AX  =  0.  We  then  set  det(A)  =  0 
to  obtain  a  relation  between  the  numerical  wavenumber  k/\  and  the  frequency  lo.  Solving  for 
k a  in  this  relation,  and  writing  in  a  form  comparable  to  (5.36),  we  have  that  the  numerical 
dispersion  relation  for  the  FEM  scheme  applied  to  Lorentz  media  is 


A:  A  —  AfeC^) 


— —  sin 

Ax: 


loa  A  2  rjr 
—  ^r\/er,F  E  > 


(5.39) 
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where  the  discrete  relative  complex  permittivity  is  given  to  be 

L  _  wi  -  es,A^o,A  +  *^awa 

er,FE  2  2  i  FT  ■ 

WA  —  W0,A  +  MAWA 


Here 


es 

^S,A  9  /99 

Aa  =  Acos(wAf/2)a 
wo, A  =  W0  cos(wAf/2)/?a2 
wa  =  ujsuql, 


(5.40) 


(5.41) 

(5.42) 

(5.43) 

(5.44) 


where  sw,  a,  (3  are  as  defined  in  equations  (5.23)— (5.25).  Thus,  the  FEM  scheme  misrepre¬ 
sents  es,  A,  wo  and  w  as  es,A,  Aa,  wo,a  and  wa  respectively.  In  particular,  note  that  for  the 
finite  element  method  applied  to  the  Lorentz  model,  Aa  does  not  depend  on  the  continuous 
model  parameter  es,  however,  wo,a  does.  In  the  Debye  model,  Aa  depended  on  f3 2  whereas 
for  the  Lorentz  model,  wo,a  depends  directly  on  /3.  However,  the  discrete  parameter  appears 
as  Wq  a  in  the  dispersion  relation,  so  the  contribution  from  j3  is  again  raised  to  the  second 
power.  This  coupled  with  the  fact  that  the  regime  of  interest  is  wA t  small,  means  that 
again  the  effect  of  es  on  the  dispersion  caused  by  the  wo.a  term  is  likely  to  be  small. 

We  compare  the  FEM  scheme  for  Lorentz  media  with  two  different  finite  difference 
schemes  which  have  been  analyzed  in  [28].  For  the  scheme  in  [18],  which  we  will  refer  to  as 
the  JHT  scheme,  the  numerical  dispersion  relation  is  given  as 


kA  —  ^jhtM  —  ^7  sin  1 


w  Az 
- —Su 


c  2  ~u/  V  oJHT 

where  the  discrete  relative  complex  permittivity  for  the  JHT  scheme  is 


,L 

er-,JHT 


W, 


—  es, Awo,A  ®^AWA 

W^  —  Wq  ^  +  zAawa 


(5.45) 


(5.46) 


Here  the  discrete  representations  of  the  continuous  model  parameters  and  the  frequency  are 
given  as 


es,A  =  <7  (5.47) 

Aa  =  A—  (5.48) 

Scj 

wo,  a  =  Woy^  cos(wA  t)  (5.49) 

wa  =  ujsu,  (5.50) 


where  sL  = 


sin(wAt) 
wA  t 


The  second  finite  difference  scheme  is  presented  in  [19]  and  will  be  refered  to  as  the  KF 
scheme.  The  numerical  dispersion  relation  for  this  scheme  is  given  by 


/v 


kf(w)  =  ^sin 


-l 


w  A  z 
-—su 


c  2 


r,KF 


(5.51) 
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where  the  discrete  relative  complex  permittivity  for  the  KF  scheme  is 


L  _  <4  -  es,Awo,A  +  *^awa 
h,KF  o  2  i  vr 

k>A  —  UJq  ^  +  zAa^a 

The  discrete  representations  of  the  continuous  model  parameters  and  the  frequency  are 
given  to  be 


(5.52) 


eSjA  =  (5.53) 

Aa  =  Acos(o;Af/2)  (5.54) 

wo, a  =  cos(u;At/2)  (5.55) 

wa  =  usu.  (5.56) 


We  now  plot  the  phase  error  as  defined  in  (5.10)  for  the  finite  element  scheme  applied 
to  Lorentz  media  and  we  compare  it  with  the  phase  errors  for  the  JHT  and  the  KF  finite 
difference  schemes.  The  phase  error  is  plotted  against  values  of  u)At  in  the  range  [0, 7r] .  We 
note  that  u>At  =  27r/lVppp,  where  Nppp  is  the  number  of  points  per  period  and  is  related  to 
the  number  of  points  per  wavelength  Nppw  via 

Nppw  =  z'-ZVppp.  (5.57) 


To  generate  the  plots  below  we  have  used  the  following  values  for  the  medium  parame¬ 
ters: 


eoo  =  1  (5.58) 

es  =  2.25  (5.59) 

t  =  1.786  x  10~16  sec  (5.60) 

loq  =  4  x  1016  rad/sec.  (5.61) 


Figures  41-43  plot  the  phase  error  $  versus  ojAt  for  the  FEM  scheme  applied  to  the 
Lorentz  model,  with  various  values  of  ho  and  u.  Figures  44-49  plot  the  same  for  the  finite 
difference  methods.  Figures  50-52  plot  the  phase  error  versus  lu  for  the  FEM  scheme 
applied  to  the  Lorentz  model,  with  various  values  of  ho  and  v.  Figures  53-58  plot  the  same 
for  the  finite  difference  methods.  Note  that  the  dispersion  for  the  finite  element  method 
reduces  as  v  goes  to  1  for  all  values  of  ho. 

Figure  59  shows  the  effect  of  varying  hT  on  the  finite  element  scheme  for  v  =  1.  Figures 
60  and  61  demonstrate  this  effect  for  the  JHT  and  KF  finite  difference  methods,  respectively. 

In  Figures  62-67  we  plot  the  real  and  imaginary  parts  of  the  relative  complex  permittiv¬ 
ity  for  Lorentz  media  for  the  continuous  equations  (exact  values),  the  finite  element  scheme 
with  v  =  {yH/2, 1,4},  ho  =  {0.1,0.01},  and  finally  the  finite  difference  scheme  with  v  =  1 
fixed  and  ho  =  {0.1,0.01}.  For  ho  =  0.1,  as  v  is  increased,  the  discrete  permittivites  of  the 
finite  element  scheme  approach  those  of  the  finite  difference  method.  However,  for  v  =  1 
it  seems  that  the  finite  element  approximation  is  actually  better.  For  ho  =  0.01,  the  agree¬ 
ment  of  the  discrete  permittivities  with  the  exact  values  is  better  than  with  ho  =  0.1,  for 
each  value  of  v.  For  v  =  1  we  see  the  best  agreement  of  the  discrete  real  and  imaginary 
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co  At  co  At 


Figure  41:  Plot  of  the  phase  error  </>  versus  ujAt  for  the  finite  element  scheme  of  the  Lorentz 
model  with  v  =  {1,  y^2/3,  ^/l/2,  .6}  using  ho  =  0.1.  (Right  plot  is  zoom  of  left  plot.) 


Lorentz  dispersion  for  FEM  with  h  =0.01 


Figure  42:  Plot  of  the  phase  error  cj)  versus  u)At  for  the  finite  element  scheme  of  the  Lorentz 
model  with  v  =  {1,  2/3,  i/l/2,  .6}  using  ho  =  0.01.  (Right  plot  is  zoom  of  left  plot.) 


permittivities  with  the  exact  values,  regardless  of  ho-  This  is  not  seen  in  the  plots  of  the 
discrete  values  of  the  parameters  A,  es  and  u>o  in  Figures  69,  68  and  70,  respectively.  In 
these  figures  we  see  that  the  discrete  parameters  in  the  finite  element  scheme  have  better 
agreement  with  the  exact  values  as  v  is  increased,  even  beyond  1.  However  the  discrete 
value  of  the  frequency  uj,  Figure  71,  seems  to  have  a  better  agreement  with  the  exact  value 
when  v  =  1.  As  the  discrete  complex  permittivity  matches  the  exact  more  closely  when 
v  =  1,  we  conclude  that  for  the  Lorentz  model,  the  value  of  v  that  correctly  represents 
ui  will  sufficiently  model  the  complex  permittivity.  For  the  finite  difference  scheme,  the 
discrete  parameters  do  not  depend  on  the  value  of  v.  However,  as  with  the  finite  element 
method,  if  ho  is  decreased  they  agree  better  with  the  exact  values. 
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co  At  co  At 


Figure  43:  Plot  of  the  phase  error  cj)  versus  uiAt  for  the  finite  element  scheme  of  the  Lorentz 
model  with  v  =  {1,  -y/ 2/3,  -^/l/2,  .6}  using  ho  =  0.001.  (Right  plot  is  zoom  of  left  plot.) 


Lorentz  dispersion  for  JHT  with  hQ=0.1  Lorentz  dispersion  for  JHT  with  hQ=0.1 


Figure  44:  Plot  of  the  phase  error  cj)  versus  ujAt  for  the  JHT  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  -^2/3,  -y/ 1/2,  .6}  using  ho  =  0.1.  (Right  plot  is  zoom  of  left 
plot.) 
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Lorentz  dispersion  for  JHT  with  h  =0.01 


Figure  45:  Plot  of  the  phase  error  cj)  versus  uiAt  for  the  JHT  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  ^2/3,  -y/l/2,  .6}  using  ho  =  0.01.  (Right  plot  is  zoom  of  left 
plot.) 


Lorentz  dispersion  for  JHT  with  hQ=0.001 


Lorentz  dispersion  for  JHT  with  hQ=0.001 


Figure  46:  Plot  of  the  phase  error  cj)  versus  ujAt  for  the  JHT  finite  difference  scheme  of  the 
Lorentz  model  with  u  =  {1,  a^/ 2/3,  ^1/2,  .6}  using  ho  =  0.001.  (Right  plot  is  zoom  of  left 
plot.) 
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Lorentz  dispersion  for  KF  with  h  =0.1 


Figure  47:  Plot  of  the  phase  error  cj)  versus  a; At  for  the  KF  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  -\/ 2/3,  -y/ 1/2,  .6}  using  ho  =  0.1.  (Right  plot  is  zoom  of  left 
plot.) 
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Lorentz  dispersion  for  KF  with  hQ=0.01 


Figure  48:  Plot  of  the  phase  error  cj)  versus  uAt  for  the  KF  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  y/2/3,  y/l/2,  .6}  using  ho  =  0.01.  (Right  plot  is  zoom  of  left 
plot.) 
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Lorentz  dispersion  for  KF  with  hQ=0.001  Lorentz  dispersion  for  KF  with  hQ=0.001 


Figure  49:  Plot  of  the  phase  error  0  versus  a; At  for  the  KF  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  y/ 2/3,  y^l/2,  .6}  using  ho  =  0.001.  (Right  plot  is  zoom  of  left 
plot.) 
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Figure  50:  Plot  of  the  phase  error  <j)  versus  ui  for  the  finite  element  scheme  of  the  Lorentz 
model  with  u  =  {1,  2/3,  i/l/2,  .6}  using  /iq  =  0.1.  (Right  plot  is  zoom  of  left  plot.) 


Figure  51:  Plot  of  the  phase  error  <j)  versus  ui  for  the  finite  element  scheme  of  the  Lorentz 
model  with  v  =  {1,  a/2/3,  a^/  1/2,  .6}  using  ho  =  0.01.  (Right  plot  is  zoom  of  left  plot.) 
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Figure  52:  Plot  of  the  phase  error  (j>  versus  ui  for  the  finite  element  scheme  of  the  Lorentz 
model  with  v  =  {1,  -y/ 2/3,  ^/l/2,  .6}  using  h$  =  0.001.  (Right  plot  is  zoom  of  left  plot.) 
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Figure  53:  Plot  of  the  phase  error  (f>  versus  to  for  the  JHT  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  -\/ 2/3,  -y/ 1/2,  .6}  using  ho  =  0.1.  (Right  plot  is  zoom  of  left 
plot.) 
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Figure  54:  Plot  of  the  phase  error  (j>  versus  c u  for  the  JHT  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  2/3,  yj  1/2,  .6}  using  ho  =  0.01.  (Right  plot  is  zoom  of  left 

plot.) 
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Figure  55:  Plot  of  the  phase  error  <j>  versus  to  for  the  JHT  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  yj 2/3,  y/l/2,  .6}  using  ho  =  0.001.  (Right  plot  is  zoom  of  left 
plot.) 
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Figure  56:  Plot  of  the  phase  error  cj)  versus  uj  for  the  KF  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  2/3,  y/l/2,  .6}  using  ho  =  0.1.  (Right  plot  is  zoom  of  left 

plot.) 


Figure  57:  Plot  of  the  phase  error  cj)  versus  uj  for  the  KF  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  y/2/3,  y/l/2,  .6}  using  ho  =  0.01.  (Right  plot  is  zoom  of  left 
plot.) 
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Lorentz  dispersion  for  KF  with  hQ=0.001 


Lorentz  dispersion  for  KF  with  h  =0.001 
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Figure  58:  Plot  of  the  phase  error  <f>  versus  uj  for  the  KF  finite  difference  scheme  of  the 
Lorentz  model  with  v  =  {1,  a^/ 2/3,  1/2,  .6}  using  ho  =  0.001.  (Right  plot  is  zoom  of  left 
plot.) 


Figure  59:  Plot  of  the  log  of  the  phase  error  cj)  versus  u)  for  the  finite  element  scheme  of  the 
Lorentz  model  with  v  =  1  using  ho  =  {0.1,  0.01,  0.001}.  (Right  plot  is  zoom  of  left  plot.) 
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Figure  60:  Plot  of  the  log  of  the  phase  error  cj)  versus  co  for  the  JHT  finite  difference  scheme 
of  the  Lorentz  model  with  v  =  1  using  ho  =  {0.1,0.01,0.001}.  (Right  plot  is  zoom  of  left 
plot.) 


Figure  61:  Plot  of  the  log  of  the  phase  error  <j)  versus  u  for  the  KF  finite  difference  scheme 
of  the  Lorentz  model  with  u  =  1  using  ho  =  {0.1,0.01,0.001}.  (Right  plot  is  zoom  of  left 
plot.) 
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Lorentz  with  h  =  0.1 ,  v  =  V(1/2) 


Figure  62:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Lorentz  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  y^l/2,  h$  = 
0.1  and  the  KF  finite  difference  scheme  with  v  =  1,  ho  =  0.1. 


Lorentz  with  hQ  =  0.01 ,  v  =  V(1/2) 


Figure  63:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Lorentz  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  y^l/2,  ho  = 
0.01  and  the  KF  finite  difference  scheme  with  v  =  1,  ho  =  0.01. 
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Lorentz  with  hQ  =  0.1 ,  v  =  1 


Figure  64:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Lorentz  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  1,  ho  =  0.1 
and  the  KF  finite  difference  scheme  with  v  =  1,  ho  =  0.1. 


Lorentz  with  h  =  0.01 ,  v  =  1 


Figure  65:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Lorentz  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  1,  ho  =  0.01 
and  the  KF  finite  difference  scheme  with  v  =  1,  ho  =  0.01. 


57 


Lorentz  with  hQ  =  0.1 ,  v  =  4 


Figure  66:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Lorentz  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  4,  ho  =  0.1 
and  the  KF  finite  difference  scheme  with  v  =  1  ,ho  =  0.1. 


Lorentz  with  h  =  0.01 ,  v  =  4 


Figure  67:  Plots  of  the  real  and  imaginary  parts  of  the  relative  complex  permittivity  for 
Lorentz  media  for  the  continuous  equations,  the  finite  element  scheme  with  v  =  4,  /io  =  0.01 
and  the  KF  finite  difference  scheme  with  v  =  1,  ho  =  0.01. 
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Lorentz  with  hQ  =  0.1  Lorentz  with  h0  =  0.01 
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Figure  68:  Plots  of  es  for  Lorentz  media  for  the  continuous  equations,  the  finite  element 
scheme  and  the  KF  finite  difference  scheme  {v  =  1)  with  ho  =  0.1  (left)  and  ho  =  0.01 
(right). 


Figure  69:  Plots  of  A  for  Lorentz  media  for  the  continuous  equations,  the  finite  element 
scheme  and  the  KF  finite  difference  scheme  {y  =  1)  with  ho  =  0.1  (left)  and  ho  =  0.01 
(right). 
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Figure  70:  Plots  of  the  resonant  frequency  ujq  for  Lorentz  media  for  the  continuous  equations, 
the  finite  element  scheme  and  the  KF  finite  difference  scheme  ( v  =  1)  with  ho  =  0.1  (left) 
and  ho  =  0.01  (right). 


Figure  71:  Plots  of  the  frequency  ui  for  Lorentz  media  for  the  continuous  equations,  the 
finite  element  scheme  and  the  KF  finite  difference  scheme  (u  =  1)  with  ho  =  0.1  (left)  and 
ho  =  0.01  (right). 
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Interrogating  signal  for  Debye  at  f=1  el  0  Hz 


Figure  72:  The  interrogating  signal  when  the  frequency  is  10  GHz. 


6  Numerical  Simulations 

To  verify  the  observations  drawn  from  the  stability  and  dispersion  analyses  in  the  previous 
sections,  we  examine  simulations  of  several  sample  problems,  which  use  either  the  Debye 
polarization  model  or  the  Lorentz  model.  In  each  case  we  will  successively  reduce  the  value 
of  hT,  and  hence  the  time  step,  holding  all  other  parameters  fixed,  until  convergence  is 
achieved. 

6.1  Debye  Simulations 

In  our  first  numerical  calculation  we  simulate  the  propagation  of  12  cycles  of  a  truncated 
sine  wave  with  carrier  frequency  at  10  GHz  (1  x  1010  Hz),  plotted  in  Figure  72,  which 
is  normally  incident  on  a  Debye  medium  from  a  vacuum.  The  medium  is  defined  by  the 
parameters  given  in  Section  4.2,  namely  es  =  78.2,  =  1,  and  r  =  8.1  x  10~12  seconds. 

We  performed  simulations  with  the  finite  element  scheme  using  v  =  4,  as  the  dispersion 
analysis  demonstrated  that  for  this  value  of  v  we  obtained  the  least  phase  error  in  the 
regime  that  is  of  interest  in  this  simulation.  A  time  trace  of  the  electric  field  at  a  depth  of 
15  mm  was  recorded.  Figure  73  shows  the  solutions  for  three  different  values  of  hT.  The 
traces  for  hT  =  0.16  and  hr  =  0.08  appear  to  be  identical  at  this  scale.  However,  Figure 
74,  which  displays  a  closer  view  of  the  central  portion  of  the  signal,  demonstrates  that  true 
convergence  has  not  been  reached.  Figure  75  displays  plots  of  one  of  the  peaks  from  the 
central  portion  of  the  signal  with  smaller  hT  values.  From  this  plot  we  can  see  that  at  most 
hT  =  0.01  (i.e. ,  at  least  100  points  per  r)  is  necessary  for  convergence  on  this  small  scale. 
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FEM  on  Debye  with  v=4,  f=1  el  0 


Figure  73:  Time  trace  of  the  electric  field  at  a  depth  of  15  mm  into  a  Debye  medium 
(frequency  is  10  GHz)  computed  using  the  finite  element  scheme,  for  three  different  values 
of  hT. 


We  note  that  while  this  frequency  is  “relatively”  small  (i.e. ,  lot  =  0.5089  <  1),  and 
therefore  the  corresponding  points  per  wavelength  is  large  in  this  example  ( ppw  =  873), 
however  this  high  resolution  of  the  wavelength  is  not  the  reason  that  the  value  of  hT  =  0.01 
is  required.  In  fact,  for  a  “relatively”  high  frequency  of  100  GHz  (i.e.,  lot  =  5.089  >  1, 
ppw  =  87.3),  we  have  essentially  the  same  restriction.  Figure  72  plots  the  interrogating 
signal  when  the  frequency  is  100  GHz.  Figure  77  shows  the  solutions  for  three  different 
values  of  hT.  Again,  the  curves  corresponding  to  the  smaller  two  values  appear  identical. 
However,  Figure  78  shows  that  at  most  hT  =  0.01  is  necessary  for  convergence  on  this 
small  scale.  In  this  example,  the  points  per  t  are  on  the  same  order  of  magnitude  as  the 
points  per  wavelength.  Any  further  increase  in  frequency  would  result  in  the  period  of  the 
interrogating  signal,  T,  being  much  less  than  the  relaxation  time,  and  therefore  chosing 
hT  =  0.01  would  not  guarantee  that  the  smallest  time  scale  in  the  problem  is  resolved. 

These  examples  verify  the  guideline  of  at  least  100  points  per  relaxation  time,  when  r 
is  the  smallest  time  scale  to  be  resolved.  This  agrees  with  FDTD  results  from  [28],  in  that 
to  ensure  minimal  dispersion  and  dissapation  error,  hT  should  be  less  than  0.01,  preferably 
O(10-3).  As  the  same  number  of  time  steps  are  required  for  both  the  finite  element  method 
and  the  finite  difference  method,  the  fact  that  the  finite  element  method  used  here  requires 
a  linear  solve  at  each  time  step  is  an  issue.  However,  both  of  the  above  simulations  for 
interrogating  a  Debye  material  using  hT  =  0.01  took  less  than  300  seconds  on  a  3.8  GHz 
speed  processor.  The  case  when  hT  =  0.16  takes  only  .3  seconds.  In  order  to  compare 
the  stability  and  phase  errors  between  the  finite  element  method  and  the  finite  difference 
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FEM  on  Debye  with  v=4,  f=1  el  0 


Figure  74:  Closer  view  of  the  center  of  the  time  trace  of  the  electric  field  at  a  depth  of  15 
mm  into  a  Debye  medium  (frequency  is  10  GHz),  for  three  smaller  values  of  hT. 


method,  we  first  determine  that  hT  =  0.04  requires  0.3  seconds  and  hT  =  0.0015  requires 
280  seconds.  Figure  79  compares  the  stability  and  phase  errors  of  the  finite  element  method 
and  the  finite  difference  method  for  the  parameter  set  that  resulted  in  a  0.3  second  runtime. 
Figure  80  displays  the  same  comparisons  for  the  280  second  runtime  case.  The  finite  element 
method  has  the  lesser  dispersion  error  only  for  the  low  frequency  and  low  resolution  case. 
However,  the  stability  results  and  phase  errors  for  the  other  cases  are  all  on  the  same  order 
of  magnitude  as  the  finite  difference  method. 
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FEM  on  Debye  with  v=4,  f=1  el  0 


Figure  75:  Magnified  view  of  the  time  trace  of  the  electric  field  at  a  depth  of  15  mm  into  a 
Debye  medium  (frequency  is  10  GHz).  Convergence  on  this  scale  is  achieved  at  hT  =  0.01. 


Interrogating  signal  for  Debye  at  f=10e10  Hz 


Figure  76:  The  interrogating  signal  when  the  frequency  is  100  GHz. 
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FEM  on  Debye  with  v=4,  f=1  Oel  0 


Figure  77:  Time  trace  of  the  electric  field  at  a  depth  of  15  mm  into  a  Debye  medium 
(frequency  is  100  GHz)  computed  using  the  finite  element  scheme,  for  three  different  values 
of  hT. 


FEM  on  Debye  with  v=4,  f=1  Oel  0 


Figure  78:  Magnified  view  of  the  time  trace  of  the  electric  field  at  a  depth  of  15  mm  into  a 
Debye  medium  (frequency  is  100  GHz).  Convergence  on  this  scale  is  achieved  at  hT  =  0.01. 
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1.05 


Debye  stability  comparison  (runtime  0.3  s) 


Figure  79:  The  stability  and  phase  errors  of  the  finite  element  method  and  the  finite  differ¬ 
ence  method  applied  to  the  Debye  model  for  the  parameter  set  that  resulted  in  a  0.3  second 
runtime. 


Debye  stability  comparison  (runtime  280  s) 


Figure  80:  The  stability  and  phase  errors  of  the  finite  element  method  and  the  finite  dif¬ 
ference  method  applied  to  the  Debye  model  for  the  parameter  set  that  resulted  in  a  280 
second  runtime. 
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Interrogating  signal  for  Lorentz  at  f=15e14  Hz 


Figure  81:  The  interrogating  signal  when  the  frequency  is  1.5  PHz. 


6.2  Lorentz  Simulations 

For  interrogation  of  a  Lorentz  medium,  we  simulate  the  propagation  of  12  cycles  of  a 
truncated  sine  wave  with  carrier  frequency  at  1.5  PHz  (1.5  x  1015  Hz),  plotted  in  Figure 
81,  which  is  normally  incident  on  the  medium  from  a  vacuum.  The  medium  is  defined  by 
the  parameters  given  in  Section  4.3,  namely  es  =  2.25,  £oo  =  1,  r  =  1.786  x  10~16,  and 
loq  =  4  x  1016.  We  performed  simulations  with  the  finite  element  scheme  using  u  =  1.  A 
time  trace  of  the  electric  field  at  a  depth  of  0.01  mm  was  recorded.  Figure  82  shows  the 
solutions  for  three  different  values  of  hT.  The  traces  for  h o  =  0.16  and  ho  =  0.08  appear  to 
be  identical  at  this  scale.  However,  Figure  83,  which  displays  a  closer  view  of  the  central 
portion  of  the  signal,  demonstrates  that  true  convergence  has  not  been  reached.  Figure  84 
displays  plots  of  one  of  the  peaks  from  the  central  portion  of  the  signal  with  smaller  ho 
values.  From  this  plot  we  can  see  that  at  most  ho  =  0.02  (i.e. ,  at  least  50  points  per  uo)  is 
necessary  for  convergence  on  this  small  scale. 

The  required  value  of  ho  <  0.02  is  in  line  with  the  suggested  range  of  0(  10~2).  This 
also  agrees  with  FDTD  results  from  [28]. 

As  we  have  simulated  propagation  deeper  (w.r.t  wavelengths)  into  the  material  in  this 
case  as  compared  with  the  Debye  simulations  above  (6600  spatial  steps,  i.e.  50  wavelengths 
into  the  material  versus  437  spatial  steps  or  0.5  wavelength  for  the  Debye  medium),  our 
computational  time  is  appropriately  longer.  For  ho  =  0.02,  the  runtime  was  greater  than 
10  minutes  on  a  3.8  GHz  speed  processor. 

Again,  the  finite  element  method  used  here  requires  a  linear  solve  at  each  time  step  as 
opposed  to  both  finite  difference  methods.  Thus,  in  order  to  compare  the  stability  plots 
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FEM  on  Lorentz  with  v=1 


Figure  82:  Time  trace  of  the  electric  field  at  a  depth  of  0.01  mm  into  a  Lorentz  medium 
computed  using  the  finite  element  scheme,  for  three  different  values  of  ho- 


and  phase  errors  of  the  finite  element  method  to  the  finite  difference  method  (here,  the 
KF  scheme),  we  first  determine  the  value  of  ho  in  the  KF  scheme  which  gives  runtimes 
comparable  to  those  of  the  cases  with  h o  =  0.16  and  ho  =  0.02  when  used  in  the  finite 
element  method.  The  runtimes  are  10  seconds  and  780  seconds,  respectively,  and  the 
corresponding  ho  values  are  0.03  and  0.0038  on  a  3.8  GHz  speed  processor.  Figure  85 
compares  the  stability  and  phase  errors  of  the  finite  element  method  and  the  finite  difference 
method  for  the  parameter  set  that  resulted  in  a  10  second  runtime.  Figure  86  displays  the 
same  comparisons  for  the  780  second  runtime  case.  Similar  to  the  Debye  case,  the  stability 
results  and  phase  errors  are  all  on  the  same  order  of  magnitude  as  in  the  finite  difference 
method. 
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FEM  on  Lorentz  with  v=1 


Figure  83:  Closer  view  of  the  center  of  the  time  trace  of  the  electric  field  at  a  depth  of  0.01 
mm  into  a  Lorentz  medium,  for  three  smaller  values  of  /iq. 


FEM  on  Lorentz  with  v=1 


Figure  84:  Magnified  view  of  the  time  trace  of  the  electric  field  at  a  depth  of  0.01  mm  into 
a  Lorentz  medium.  Convergence  on  this  scale  is  achieved  at  hT  =  0.02. 
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1.001 


Lorentz  stability  comparison  (runtime  10  s) 


Figure  85:  The  stability  and  phase  errors  of  the  finite  element  method  and  the  finite  dif¬ 
ference  method  applied  to  the  Lorentz  model  for  the  parameter  set  that  resulted  in  a  10 
second  runtime. 


Lorentz  stability  comparison  (runtime  780  s) 
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Figure  86:  The  stability  and  phase  errors  of  the  finite  element  method  and  the  finite  dif¬ 
ference  method  applied  to  the  Lorentz  model  for  the  parameter  set  that  resulted  in  a  780 
second  runtime. 
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7  Conclusions 


From  the  stability  and  dispersion  analysis  as  well  as  the  simulations  shown  in  the  paper,  we 
can  conclude  that  the  artificial  dissipation  in  the  finite  element  schemes  presented  here  for 
Debye  and  Lorentz  media  are  strongly  dependent  on  the  quantity  hT  when  r  is  the  smallest 
time  scale.  For  Lorentz  media,  in  addition,  the  quantity  ho  may  be  the  dominant  quantity 
if  To  =  27r/o;o  is  smaller  than  r.  We  see  that  hT  or  ho  have  to  be  sufficiently  small  so  that 
not  much  energy  is  lost  after  the  large  amount  of  timesteps  usually  needed  to  propagate 
short  pulses  any  appreciable  distance  into  a  medium.  For  Debye  media  hT  is  recommended 
to  be  at  least  100  points  per  r,  preferably  hT  =  O(10~3).  For  Lorentz  media  we  recommend 
either  ho  or  hT  to  be  0(  10-2)  in  order  to  minimize  dissipation.  These  results  agree  with 
the  guidelines  posed  in  [28]  for  FDTD  schemes. 

From  the  dispersion  analysis  for  the  finite  element  schemes  we  see  that  the  value  of  v 
to  be  chosen  should  be  the  one  that  provides  the  best  agreement  of  the  discrete  relative 
complex  permittivity  with  the  exact  complex  permittivity,  i.e.  the  value  that  has  the  least 
phase  error  for  the  regime  of  interest.  The  unconditional  stability  of  the  finite  element 
scheme  allows  the  user  to  have  some  freedom  in  the  choice  of  v.  For  the  Debye  model,  the 
value  of  v  that  best  represents  A  should  sufficiently  model  the  complex  permittivity.  For 
the  Lorentz  model,  the  value  of  v  that  correctly  represents  lo  should  be  chosen. 
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